一種基于兩階段投影調(diào)整的快速ct圖像重建方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種基于兩階段投影調(diào)整的快速CT圖像重建方法，屬于生物醫(yī)學(xué)成 像、無損檢測等技術(shù)領(lǐng)域。
【背景技術(shù)】
[0002] 計(jì)算機(jī)斷層掃描（CT)成像是一種重要的無損檢測技術(shù)，經(jīng)過幾十年的發(fā)展，已經(jīng) 從醫(yī)學(xué)領(lǐng)域擴(kuò)展到工業(yè)等其他領(lǐng)域。隨著臨床應(yīng)用與工業(yè)應(yīng)用范圍的日益擴(kuò)大與深入，受 檢測環(huán)境、成本、時(shí)間、人員健康等諸多因素限制，不完全投影數(shù)據(jù)下的CT圖像重建已不可 避免。在現(xiàn)今主流商用CT機(jī)中，以濾波反投影算法(FBP)為代表的解析法是最常用的重建算 法。該類算法的優(yōu)點(diǎn)是計(jì)算量小，重建速度快，但要求獲得充足的投影數(shù)據(jù)。在不完全投影 情況下，現(xiàn)有的濾波反投影類算法難以重建出高質(zhì)量的CT圖像。
[0003] 與解析法不同，代數(shù)重建算法將CT圖像重建問題轉(zhuǎn)化為大規(guī)模線性方程組求解問 題，通過結(jié)合一些正則化先驗(yàn)知識，在投影數(shù)據(jù)不完全時(shí)仍可重建高質(zhì)量圖像。另外，代數(shù) 重建算法不易受所處理問題模型的限制并且可以引入被檢物體的先驗(yàn)信息，容易擴(kuò)展到CT 應(yīng)用的其他領(lǐng)域，展示出了很大的發(fā)展空間與應(yīng)用潛力。
[0004] 雖然世界上第一臺CT機(jī)所采用的圖像重建算法就是代數(shù)重建法，但計(jì)算量大、重 建速度慢的缺點(diǎn)，在很長一段時(shí)間內(nèi)限制了代數(shù)重建法在CT成像系統(tǒng)的大規(guī)模應(yīng)用。目前 隨著并行計(jì)算理論及計(jì)算機(jī)技術(shù)的快速發(fā)展，代數(shù)重建法的技術(shù)瓶頸已得到很大程度的緩 解，其優(yōu)勢更加凸顯，近年來代數(shù)重建技術(shù)重新引起了眾多研究者的注意。如何在不完全投 影情況下提高代數(shù)重建法的收斂速度與精度，成為了代數(shù)重建算法能否在實(shí)際成像系統(tǒng)中 得到普遍應(yīng)用的關(guān)鍵。

【發(fā)明內(nèi)容】

[0005] 本發(fā)明提供一種基于兩階段投影調(diào)整的快速CT圖像重建方法，達(dá)到收斂速度快， 能夠精確求解非一致線性方程組的最小二乘解，可用于不完全投影情況下的快速CT圖像重 建的目的。
[0006] 本發(fā)明的理論分析為：
[0007] 代數(shù)重建法(ART)將圖像重建問題轉(zhuǎn)化為如下大規(guī)模線性方程組求解問題：
[0008] Ax = b (1)
[0009] 其中，X即為待求的圖像，為NX 1向量，由對二維或三維圖像按一維重新排列得到， X的每個(gè)元素 Xl對應(yīng)原圖像的一個(gè)像素。A是Μ X N系統(tǒng)矩陣，每個(gè)元素描述第j個(gè)像素對第 i條射線的貢獻(xiàn)。b為MX 1測量數(shù)據(jù)向量，其第i個(gè)元素 h對應(yīng)第i條射線的衰減。
[00?0] 最早的ART方法由Kaczmarz提出，也稱Kaczmarz方法（簡稱KA)。其基本思想是：將 各方程看作N維空間的一個(gè)超平面，從初始解向量X*3出發(fā)，依次在各超平面上投影，逐步逼 近方程組的真實(shí)解，可用如下數(shù)學(xué)公式表示：
[0011] xk+1 = P(b,xk) =PM〇---〇P2°Pi(b,xk) (2)
[0013] .(4)
[0012] (3)
[0014] (5)
[0015] 其中，ai為矩陣A的第i行向量，〈·〉與Μ · I I分別表示內(nèi)積與歐拉范數(shù)，aiT表示對 的轉(zhuǎn)置表示將向量投影到第i個(gè)超平面的操作，表示KA第k次迭代時(shí)投 影到第j個(gè)超平面的向量。每次迭代，對向量xki行操作P(b，xk)，即依次在各超平面上進(jìn)行 一輪投影，得到下一代的解向氣
[0019] (8)[0020] 式(7)中的參數(shù)β根據(jù)如下表達(dá)式計(jì)算：
[0016] 為了提高Kaczmarz方法的效率，我們提出了基于投影調(diào)整的快速代數(shù)重建算法。 在每次迭代對KA投影到超平面i的投影向量x 1^1沿其與前代該超平面投影向量，-1，1連線方 向進(jìn)行進(jìn)一步的調(diào)整優(yōu)化，使得調(diào)整后的新向量，在一― 1與一"連線方向上最優(yōu)(即到真實(shí) 解向量f的平方距離或誤差最?。；谕队罢{(diào)整的快速代數(shù)重建算法每次迭代的數(shù)學(xué)公 式如下：
[0017] (：6)
[0018] (7)
[0021] (9)
[0022]
[0023] (11)
[0024]其中，d15，1表示向量X15，1到真實(shí)解向量f的平方誤差（或距離）。在ΚΑ中，將向量x k，j 投影到第i個(gè)平面得到χΜ，對應(yīng)的必1滿足如下關(guān)系：
[0025]
(12)
[0026] 這樣由一個(gè)初始值f開始((^二#)，根據(jù)式（12)依次迭代計(jì)算d1-1，1與必 1，再將 dk~與必1代入式（10)，可消去d°，求得Ad。將式（10)計(jì)算得到的Ad代入式(9)求出β，再將 求出的β代入式⑴即可求出調(diào)整后的Α然后根據(jù)式(13)、（14)更新與d M:
[0027] (13)
[0028] (14)
[0029] 上述快速代數(shù)重建算法中每次迭代可以采用不同的超平面投影順序。在每次迭代 采用隨機(jī)的投影順序，許多情況下可以明顯提高算法的效率。本發(fā)明后面的例子中，也將采 用隨機(jī)投影順序的快速代數(shù)重建算法。
[0030] 理論分析表明，如果待求解線性方程組（1)是一致線性方程組（即方程組存在真實(shí) 解），則KA將收斂到方程組的解。若初始解/屬于Α τ的值域R(AT)，KA收斂的解是方程組的最 小二范數(shù)解。若方程組（1)是非一致線性方程組（即方程組不存在解），則KA最終只能收斂到 一個(gè)區(qū)域，而不是某個(gè)特定解?；谕队罢{(diào)整的快速代數(shù)重建法在求解一致線性方程組的 收斂性與KA相同，在非一致方程組情況下也不能收斂到某個(gè)特定解。
[0031]由于噪聲及正則化的影響，在很多情況下基于代數(shù)重建法的CT圖像重建問題被轉(zhuǎn) 化為非一致線性方程組求解問題。由于非一致線性方程組Ax = b無解，我們轉(zhuǎn)而求該方程組 均方誤差I(lǐng) |Ax-b| |2最小的解，即最小二乘解。為此，本發(fā)明將非一致線性方程組的求解轉(zhuǎn) 化為兩個(gè)一致線性方程組的求解，并利用基于投影調(diào)整的快速代數(shù)重建算法分別求解兩個(gè) 一致線性方程組，最終實(shí)現(xiàn)非一致線性方程組最小二乘解的精確求解與CT圖像的快速重 建。其理論依據(jù)與實(shí)施方法如下。
[0032]若方程組（1)為非一致線性方程組，則將方程組右邊的向量b分解為互相垂直的兩 個(gè)分量：
[0033] b = bi<A)+bAr(A^ ( 15)
[0034] 其中，bR(A)為b在矩陣A的列空間或值域空間R(A)的投影分量，bv(A/) Sb在矩陣Ατ 的零空間N(AT)的投影分量。假設(shè)/為非一致線性方程組的最小二乘解，則f滿足：
[0035] Ax* = bR(A) (16)
[0036] 而實(shí)際是如下方程的解：
[0037] ATy = 〇 (17)
[0038] 其中，y為MX 1向量。方程組（16)與（17)都是一致線性方程組，都可以用基于投影 調(diào)整的快速代數(shù)重建算法求精確解。因此我們可以以b作為y的初始解y* 3,先求出式（17)的 解\^，然后利用bw=b-求出bR(A)，最后再解方程組（16)求出最小二乘解X'為 加速求解過程，我們交互地迭代求解y與X。每次迭代的數(shù)學(xué)公式如下：
[0039] (18)
[0040] (19)
[0041] (20)
[0042] 在算法起始收斂階段，由式(20)迭代產(chǎn)生的xk+1實(shí)際上是收斂到方程組(21)的解 而不是方程組（16)對應(yīng)的最小二乘解，。
[0043] Ax = bk+1 (21)
[0044] 但是，隨著yk+1逐漸收斂到b#(Ar)，bk+1將逐漸收斂到b R(A)。這樣方程組(21)就變?yōu)?了（16)，因此xk最終仍可以收斂到與（16)對應(yīng)的最小二乘解X'但是在算法起始收斂階段， xk的收斂方向可能偏離X'
[0045] 為進(jìn)一步減小算法計(jì)算量與加快算法收斂速度，本發(fā)明將算法