【正文】 重建圖像與原圖相差很遠(yuǎn)，不能達(dá)到我們預(yù)期的結(jié)果，也不能有效地抑制噪聲。所以必須要將圖像近似地利用有限個(gè)參數(shù)來表示，這樣才能利用有限的已知來求解這有限個(gè)未知的參數(shù)來得到重建圖像。我們可以利用級(jí)數(shù)展開的數(shù)學(xué)方法將連續(xù)的近似表示為有限個(gè)參數(shù)和基函數(shù) (其中為基的個(gè)數(shù)，為索引值)的加權(quán)求和的形式。通過式()的變換可以看出已經(jīng)將連續(xù)的無限模型變成了有限的離散模型。最常用的是利用以像素為單位的基函數(shù)來表示一幅二維圖像。物理模型將直接影響到重建過程的計(jì)算速度和圖像的質(zhì)量。為此應(yīng)該在圖像質(zhì)量和計(jì)算的速度之間進(jìn)行均衡的選取物理模型。從投影和反投影的角度來看，系統(tǒng)矩陣正好反映了投影值和成像截面之間的關(guān)系，也就是成像系統(tǒng)的模型，并且是表示投影方式的函數(shù)。 系統(tǒng)矩陣的求解可先將整個(gè)圖像區(qū)域劃分為的有限個(gè)像素，將個(gè)正方形網(wǎng)格覆蓋在圖像上，設(shè)網(wǎng)格長(zhǎng)和寬皆為1。將射線看成寬為，間隔為的平行射束。系統(tǒng)矩陣中的各個(gè)(表示第個(gè)像素對(duì)第條射線投影的貢獻(xiàn))的取值遵循如下規(guī)則：，號(hào)射線通過號(hào)像素內(nèi)任一點(diǎn) ，其他 ()n*nn…21 射線與網(wǎng)格相交關(guān)系示意圖，以圖像中心為坐標(biāo)原點(diǎn)，水平和豎直方向?yàn)樽鴺?biāo)軸建立二維直角坐標(biāo)系。為方便討論，考慮射線斜率的情況，對(duì)于和的情況，可利用對(duì)稱的性質(zhì)得出。 ()然后計(jì)算與射線的相交像素的編號(hào)；將射線與垂直網(wǎng)格線交點(diǎn)橫坐標(biāo)，與水平網(wǎng)格線交點(diǎn)橫坐標(biāo)，按從小到大的順序合并排列，以數(shù)組形式保存，設(shè)為數(shù)組。顯然射線與網(wǎng)格線最多有個(gè)交點(diǎn)，即數(shù)組、只需開辟的空間。由以上計(jì)算，分四種情況討論：(1) 當(dāng)時(shí)，與射線相交網(wǎng)格編號(hào)的表達(dá)式為：； ()(2) 當(dāng)時(shí)，與射線相交網(wǎng)格編號(hào)的表達(dá)式為：,其中表示向下取整； ()(3) 當(dāng)時(shí)，與射線相交網(wǎng)格編號(hào)的表達(dá)式為：； ()(4) 當(dāng)時(shí)，與射線相交網(wǎng)格編號(hào)的表達(dá)式為：。編號(hào)以數(shù)組形式保存。只需開辟一個(gè)空間數(shù)組，設(shè)為數(shù)組[14]。但目前ART尚存在計(jì)算量大，重建速度慢，無法用硬件實(shí)現(xiàn)等缺點(diǎn)，限制了其進(jìn)一步應(yīng)用[16]。先假設(shè)一初始圖像，然后根據(jù)求一次近似圖像，再根據(jù)求二次近似圖像，如此繼續(xù)，直到滿足預(yù)定條件而后止。只考慮一條射線(例如號(hào)射線)的射線投影響，所修正的像素值也限于號(hào)射線經(jīng)過的那些像素?？傊看涡Ｕ来慰紤]一條射線的射線和，并校正該射線經(jīng)過的像素。但也有某些形式的ART并不對(duì)應(yīng)哪個(gè)最優(yōu)準(zhǔn)則，然而行之有效，本論文中討論ARTⅡ就屬于這一類。 ()利用線性不等式方程組的松弛法可得其解為：， ，任意； ，其他 () ()上式中代表迭代次數(shù)，代表取整數(shù)。 ARTⅡ算法的計(jì)算機(jī)實(shí)現(xiàn)步驟利用ARTⅡ算法重建過程如下：(1) 確定物理模型，計(jì)算整個(gè)系統(tǒng)的投影矩陣；(2) 設(shè)定迭代次數(shù)，初始化圖像矢量，對(duì)每一個(gè)圖像象素利用式()進(jìn)行校正，完成一次迭代過程；(3) 重復(fù)步驟(2)進(jìn)行下一次的迭代過程，對(duì)圖像矢量再次校正；(4) 判斷是否達(dá)到設(shè)定的迭代次數(shù)，若達(dá)到則停止迭代，否則進(jìn)入步驟2，進(jìn)行下一次迭代。本論文只研究64條平行射束條件下，18個(gè)投影角度和36個(gè)投影角度兩種情況的圖像重建結(jié)果。當(dāng)投影角度數(shù)量為18時(shí)，原始圖像、(a)(b)所示，(a)、(b)、(c)、(d) (e)(f)所示。當(dāng)投影角度數(shù)量為36時(shí)，原始圖像、(a) (b)所示，(a)、(b)、(c)、(d) (e)(f)所示。(a) 原始圖像 (b) 64*18無噪聲投影數(shù)據(jù) 原始圖像和投影數(shù)據(jù)(a)迭代20000次圖像 (b)迭代15000次圖像 (c)迭代10000次圖像(d)迭代5000次圖像 (e)迭代2000次圖像 (f)迭代1500次圖像 無噪聲時(shí)，平行射束數(shù)量為64條，18個(gè)投影角度下各重建圖像(a) 原始圖像 (b) 64*36無噪聲投影數(shù)據(jù) 原始圖像和投影數(shù)據(jù)(a) 迭代20000次圖像 (b) 迭代15000次圖像 (c) 迭代10000次圖像(d) 迭代5000次圖像 (e) 迭代2000次圖像 (f) 迭代1500次圖像 無噪聲時(shí)，平行射束數(shù)量為64條，36個(gè)投影角度下各重建圖像 有噪聲影響本文加入的噪聲同樣是以投影數(shù)據(jù)矩陣為參數(shù)的泊松噪聲。利用ARTⅡ重建算法分別迭代20000次、15000次、10000次、5000次、(a)、(b)、(c)、(d) (e)(f)所示。當(dāng)投影角度數(shù)量為36時(shí)，(a)(b)所示，(a)、(b)、(c)、(d) (e)(f)所示。(a) 原始圖像的投影數(shù)據(jù) (b) 64*18有噪聲平行投影 無噪聲時(shí)的投影數(shù)據(jù)和有噪聲影響的投影數(shù)據(jù) (a) 迭代20000次圖像 (b) 迭代15000次圖像 (c) 迭代10000次圖像(d) 迭代5000次圖像 (e) 迭代2000次圖像 (f) 迭代1500次圖像 有噪聲時(shí)，平行射束數(shù)量為64條，18個(gè)投影角度下各重建圖像由以上的重建結(jié)果可以看出，與濾波反投影算法不同，即使在射束數(shù)量和投影角度較少的條件下，ARTⅡ算法仍可以有效地進(jìn)行圖像重建，重建出原圖像的大概面貌,且有一定地抗噪聲能力，而其重建時(shí)間大于濾波反投影重建算法所用的時(shí)間。同時(shí)我們可以看到，隨著投影角度的增加，用于重建的數(shù)據(jù)量增大，由此可以得到，平行射束數(shù)量一定時(shí)，圖像重建效果隨著投影角度增加會(huì)越來越好。由于該方法對(duì)數(shù)據(jù)完備性和噪聲的要求限制了在數(shù)據(jù)不完備時(shí)的應(yīng)用，掃描速度的進(jìn)一步提高和患者所受照射劑量的降低。統(tǒng)計(jì)重建算法對(duì)噪聲的魯棒性和對(duì)數(shù)據(jù)的完備要求更優(yōu)，因此成為了現(xiàn)階段重建方法的研究方向。需要進(jìn)一步的優(yōu)化和簡(jiǎn)化算法，以提高重建的速度，而目前的研究重點(diǎn)也集中于這方面。缺點(diǎn)：(1) 計(jì)算量很大，不能實(shí)時(shí)重建；(2) 模型與解析重建和代數(shù)重建相比非常復(fù)雜；(3) 與重建模型相對(duì)應(yīng)的算法實(shí)現(xiàn)也相對(duì)復(fù)雜。已知條件為系統(tǒng)的物理模型，統(tǒng)計(jì)模型以及系統(tǒng)的輸入，輸出，需要確定的為模型中的參數(shù)。然后就是利用已有的關(guān)系確定行之有效的求解方法。在此將重點(diǎn)討論后的三個(gè)步驟。而統(tǒng)計(jì)模型的建立是在對(duì)特定成像系統(tǒng)的特性的統(tǒng)計(jì)結(jié)果和一定的先驗(yàn)知識(shí)的基礎(chǔ)之上的。本論文不對(duì)詳細(xì)的統(tǒng)計(jì)學(xué)過濾的推導(dǎo)過程做過多闡述，但是值得一提的是對(duì)于CT，這些推倒的前提都是基于X射線是服從泊松分布的。 優(yōu)化準(zhǔn)則 統(tǒng)計(jì)重建問題的“病態(tài)性”使只能在一定的優(yōu)化準(zhǔn)則下考慮問題的求解。ML準(zhǔn)則僅能在一定迭代次數(shù)范圍內(nèi)起到抵抗噪聲的作用，然而隨著迭代深入，該準(zhǔn)則不能有效的抑制噪聲，圖像質(zhì)量將變差。下面本文將對(duì)這兩個(gè)優(yōu)化準(zhǔn)則分別進(jìn)行闡釋。1982年，Shepp和Vardi將EM算法引入到PET圖像重建當(dāng)中。另一方面，EM算法是從統(tǒng)計(jì)角度入手，是一種參數(shù)估計(jì)方法，因而以概率的形式來討論由此重建的圖像的誤差問題。因此,在使用最大似然估計(jì)方法時(shí),被估計(jì)的參數(shù)假定是常數(shù)，但未知。未知參數(shù)的最大似然估計(jì)可記作: ()式中，表示隨機(jī)變量的個(gè)觀測(cè)值，令矢量形式，隨機(jī)變量的概率密度函數(shù)有未知參數(shù)決定，記為。最大似然最大期望值法圖像重建（Maximization Likelihood Expectation Maximization，MLEM）旨在尋找與測(cè)量的投影數(shù)據(jù)具有最大似然（ML）的估計(jì)解，其迭代過程是由最大期望值算法（EM）來實(shí)現(xiàn)的。這里完備數(shù)據(jù)空間定義為其中表示從第個(gè)單元發(fā)射而被第個(gè)探測(cè)器探測(cè)到的光子數(shù)。最大化的過程要比直接最大化似然函數(shù)容易。對(duì)于本文研究的圖像重建，可先將圖像區(qū)域離散化為像素的數(shù)字圖像。假設(shè)對(duì)物體斷面進(jìn)行掃描后一共獲得了個(gè)投影, 表示沿著第條投影線所探測(cè)到的投影值。這里有,。然后我們通過如下的方法來實(shí)現(xiàn)圖像的重建：給定一組測(cè)量得到的投影矢量，找到一組圖像矢量使得概率有最大值。這個(gè)Poisson分布由可以表示為： ()由于測(cè)量得到投影值只依賴于基于圖像矢量的它的期望值，因此又有： ()ML方法就是在限制為非負(fù)條件下使得這個(gè)概率最大，也就是使其似然函數(shù)大。ML類算法的迭代修正準(zhǔn)則是投影的估計(jì)值與原始采集的投影數(shù)據(jù)越接近越好，而投影的估計(jì)值是無約束的，因而會(huì)出現(xiàn)重建圖像中的噪聲不斷增加，而投影的估計(jì)值卻與原始投影數(shù)據(jù)越來越接近的情況。最大后驗(yàn)概率重建算法(maximum a posteriori, MAP)根據(jù)圖像的先驗(yàn)知識(shí)分布則可以有效的抑制噪聲。MAP重建算法是一種貝葉斯重建(Bayesian Reconstruction, BR)算法,該類算法既考慮了投影數(shù)據(jù)的統(tǒng)計(jì)特性,又充分利用了源、射線的先驗(yàn)知識(shí)分布，多次迭代后仍可以很好的抑制噪聲而且可以克服ML重建算法收斂慢的缺點(diǎn)。同MLEM算法，令矢量和矢量分別代表投影數(shù)據(jù)和重建圖像，MAP類算法采用的重建準(zhǔn)則為求在已知條件下尋求使下面的后驗(yàn)概率最大 () 其中為先驗(yàn)知識(shí)分布,以Gibbs先驗(yàn)知識(shí)為例。MAP類算法尋求式()的后驗(yàn)概率最大，等同于下式中尋求對(duì)數(shù)后驗(yàn)概率最大： ()其中為式()中分母和式()中的貢獻(xiàn)，為常數(shù)。上式和MLEM迭代公式相似，在分母上多一項(xiàng)。我們采用馬爾可夫隨機(jī)場(chǎng)(Markov random fields, MRF)描述圖像,不同的先驗(yàn)知識(shí)分布模型里又可選擇不同的勢(shì)函數(shù)。 高斯先驗(yàn)知識(shí)分布 以下是先驗(yàn)知識(shí)分布里不同的鄰域的選擇方法廣義高斯先驗(yàn)分布在鄰域類型確定后選擇不同的勢(shì)函數(shù)=零階鄰域 泊松先驗(yàn)知識(shí)分布一階鄰域吉布斯先驗(yàn)知識(shí)分布先驗(yàn)知識(shí)分布二階鄰域伽馬先驗(yàn)知識(shí)分布混合伽馬先驗(yàn)知識(shí)分布最大熵先驗(yàn)知識(shí)分布小波模型先驗(yàn)知識(shí)分布 先驗(yàn)知識(shí)分布 部分常用勢(shì)函數(shù) 抑制噪聲并能保護(hù)圖像邊緣，凹函數(shù)勢(shì)函數(shù) 特點(diǎn) 只能在抑制噪聲和保護(hù)圖像邊緣中折衷選取，凹函數(shù) 抑制噪聲并能保護(hù)圖像邊緣，凸函數(shù) 抑制噪聲并能保護(hù)圖像邊緣，能有效處理低劑量下的投射成像 ，凸函數(shù) 在MAP重建算法中采用比較多的是Gibbs先驗(yàn)知識(shí)分布。在實(shí)際重建過程中會(huì)存在如下的事實(shí)：同一幅圖象,選擇不同的Gibbs prior,則取得最佳效果的值不同；Gibbs prior函數(shù)相同,但圖像不同,則取得最佳效果的值不同；圖像相同,Gibbs prior函數(shù)相同,信噪比不同,則取得最佳效果的值不同；圖像相同，Gibbs prior函數(shù)相同，信噪比相同，但如果噪聲模型不同，則取得最佳效果的值不同[19]。 MAP算法的計(jì)算機(jī)實(shí)現(xiàn) 利用MAP算法重建過程如下：(1) 確定物理模型，計(jì)算整個(gè)系統(tǒng)的投影矩陣；(2) 設(shè)定迭代次數(shù)，初始化圖像矢量，對(duì)每一個(gè)圖像象素利用式()進(jìn)行校正，完成一次迭代過程；(3) 重復(fù)步驟(2)進(jìn)行下一次的迭代過程，對(duì)圖像矢量再次校正；(4) 判斷是否達(dá)到設(shè)定的迭代次數(shù)，若達(dá)到則停止迭代，否則進(jìn)入步驟2，進(jìn)行下一次迭代。同樣，我們將在64條平行射束條件下，18個(gè)投影角度和36個(gè)投影角度兩種情況的圖像重建結(jié)果。 無噪聲影響對(duì)于MLEM算法：設(shè)平行射束數(shù)量為64條，迭代次數(shù)為10次。當(dāng)投影角度數(shù)量為36時(shí)，(b)所示,。對(duì)于MLEM算法：設(shè)平行射束數(shù)量為64條，迭代次數(shù)為10次。當(dāng)投影角度數(shù)量為36時(shí)，(b)所示,。在像素一階鄰域內(nèi)選取勢(shì)函數(shù)，則圖像的能量函數(shù),也即懲罰函數(shù)，參數(shù)在水平方向上取值為1，Gibbs參數(shù)取不同的值[20,21]。當(dāng)投影角度數(shù)量為18時(shí)，MAP重建算法=、=、==(a)、(b)、(c)(d)所示。當(dāng)投影角度數(shù)量為36時(shí)。、。也不難看到MAP算法比MLEM算法在抗噪聲能力方面更加突出，圖像重建效果更好。 對(duì)于MAP算法而言，取不同的值，對(duì)應(yīng)得到的重建圖像效果也不同，隨著的增大，圖像中心橢圓形圖案的顏色越來越明亮，噪聲顆粒現(xiàn)象越來越模糊，圖像對(duì)比度也越強(qiáng)。 統(tǒng)計(jì)重建方法的改進(jìn)前面討論的MLEM算法存在迭代法的運(yùn)算量巨大，運(yùn)算時(shí)間過長(zhǎng)的缺點(diǎn)，這極大地限制了它的應(yīng)用。有序子集最大期望值方法（Ordered Subset