張蕾+徐伯慶

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

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

DOIDOI：10.11907/rjdk.171337

中圖分類號(hào)：TP317.4

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

0 引言

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

1 ART算法原理

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

在此，采用平行束投影方式進(jìn)行數(shù)據(jù)采集。平行射線的總數(shù)為M，圖像像素的總數(shù)為N，待重建區(qū)域由N個(gè)邊長(zhǎng)為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條射線對(duì)第j個(gè)像素的奉獻(xiàn)值，pi（1≤i≤M）為待重建圖像在第i條射線上的投影值，即待重建圖像沿此射線方向的線積分。由于M、N的值通常都很大，且每條射線與像素的交點(diǎn)個(gè)數(shù)較少，權(quán)重因子矩陣成為大型稀疏矩陣，所以只能采用迭代方式求解。

ART 算法通過(guò)對(duì)未知圖像矩陣賦初值，再將射線掃描圖像得到的實(shí)際投影數(shù)據(jù)與估計(jì)的投影數(shù)據(jù)進(jìn)行比較，使用迭代公式對(duì)初值進(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ù)訪問(wèn)順序改進(jìn)

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果與分析

3.1 圖像質(zhì)量評(píng)價(jià)方式

圖像質(zhì)量的好壞可用主觀和客觀兩種方式評(píng)價(jià)。主觀評(píng)價(jià)就是通過(guò)人眼進(jìn)行直觀判斷，這種方式雖然直觀可見，但容易受到人的主觀影響，缺乏科學(xué)依據(jù)?？陀^評(píng)價(jià)就是通過(guò)計(jì)算機(jī)模擬計(jì)算出圖像的各項(xiàng)衡量指標(biāo)，結(jié)合原始圖像和重建圖像的匹配程度來(lái)評(píng)價(jià)。為了直觀對(duì)比，本文綜合這兩種方式來(lái)評(píng)估圖像質(zhì)量，并采用以下3項(xiàng)指標(biāo)作為客觀評(píng)價(jià)標(biāo)準(zhǔn)：（1）密度曲線對(duì)比圖，即根據(jù)圖像的某一列或某一行的像素值繪制密度曲線圖，對(duì)比分析重建圖像與原始圖像像素的偏離大小，直觀判斷重建算法的優(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）歸一化平均絕對(duì)距離r，計(jì)算公式如下： r = ∑0≤i < Nxi -xi ∑0≤i < Nxi （5）由式（5）可知，r值越小，圖像失真越小，重建效果越好。

3.2 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

4 結(jié)語(yǔ)

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

參考文獻(xiàn)：

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

[2] 徐國(guó)良，陳沖，李明.圖像重構(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ù)計(jì)算的快速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é)與探測(cè)技術(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