改進(jìn)投影數(shù)據(jù)訪問順序的CT圖像重建研究

張蕾+徐伯慶

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

關(guān)鍵詞：CT圖像重建；代數(shù)重建；投影數(shù)據(jù)；訪問順序；最小相關(guān)性

DOIDOI：10.11907/rjdk.171337

中圖分類號：TP317.4

文獻(xiàn)標(biāo)識碼：A 文章編號文章編號：1672-7800（2017）008-0189-04

0 引言

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

1 ART算法原理

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

在此，采用平行束投影方式進(jìn)行數(shù)據(jù)采集。平行射線的總數(shù)為M，圖像像素的總數(shù)為N，待重建區(qū)域由N個邊長為1的正方形網(wǎng)格組成。按照CT圖像重建原理以及構(gòu)造的數(shù)學(xué)模型，投影數(shù)據(jù)與像素值之間的關(guān)系可用下式表示：

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

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

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

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

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

2 投影數(shù)據(jù)訪問順序改進(jìn)

2.1 最小相關(guān)性原理

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

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

2.2 改進(jìn)的投影數(shù)據(jù)訪問順序

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

3 實驗結(jié)果與分析

3.1 圖像質(zhì)量評價方式

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

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

3.2 實驗結(jié)果與分析

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

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

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

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

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

4 結(jié)語

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

參考文獻(xiàn)：

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

[2] 徐國良，陳沖，李明.圖像重構(gòu)的數(shù)值方法[M].北京：科學(xué)出版社，2015.

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

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

[5] 楊文良，魏東波.一種改進(jìn)投影系數(shù)計算的快速ART算法[J].CT理論與應(yīng)用研究，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] 李志鵬，叢鵬，鄔海峰.代數(shù)迭代算法進(jìn)行CT圖像重建的研究[J].核電子學(xué)與探測技術(shù)，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] 王宏鈞，路宏年，傅建.代數(shù)重建技術(shù)中投影序列選擇次序的研究[J].光學(xué)技術(shù)，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