何　騫，黃立宏1.湖南城市學(xué)院 信息科學(xué)與工程學(xué)院，湖南 益陽(yáng) 41300.湖南大學(xué) 數(shù)學(xué)與計(jì)量經(jīng)濟(jì)學(xué)院，長(zhǎng)沙 41008

結(jié)合各向異性中值擴(kuò)散的PET圖像重建算法*

何騫1,2+，黃立宏2
1.湖南城市學(xué)院 信息科學(xué)與工程學(xué)院，湖南 益陽(yáng) 413002
2.湖南大學(xué) 數(shù)學(xué)與計(jì)量經(jīng)濟(jì)學(xué)院，長(zhǎng)沙 410082

HE Qian,HUANG Lihong.PET image reconstruction algorithm combined with anisotropic median-diffusion. Journal of Frontiers of Computer Science and Technology,2016,10(8)：1166-1175.

摘要：為了有效提高正電子發(fā)射斷層成像圖像的質(zhì)量，通過(guò)把各向異性中值擴(kuò)散濾波器融合到中值根先驗(yàn)算法中，提出了一種新的基于Bayesian理論的圖像重建算法。新算法的每次迭代過(guò)程都可以分為兩步：首先用各向異性中值擴(kuò)散濾波器抑制重建圖像中的噪聲；然后用中值根先驗(yàn)算法重建圖像。仿真實(shí)驗(yàn)結(jié)果表明，在正電子發(fā)射斷層成像中，新算法不僅能有效地抑制噪聲，還能精確地保護(hù)圖像的邊緣。此外，與其他類(lèi)似算法相比，新算法吸收了各向異性中值擴(kuò)散濾波器的優(yōu)點(diǎn)，在迭代過(guò)程中對(duì)梯度閾值和擴(kuò)散次數(shù)不敏感，易于實(shí)現(xiàn)，實(shí)用性強(qiáng)。

關(guān)鍵詞：正電子發(fā)射斷層成像（PET）；各向異性中值擴(kuò)散；中值根先驗(yàn)；圖像重建；抑制噪聲

1　引言

正電子發(fā)射斷層成像（positron emission tomography，PET）是核醫(yī)學(xué)領(lǐng)域比較先進(jìn)的臨床檢查影像技術(shù)，是目前唯一可在活體上顯示生物分子代謝、受體及神經(jīng)介質(zhì)活動(dòng)的新型影像技術(shù)，現(xiàn)已廣泛地應(yīng)用于多種疾病的診斷以及新藥物的開(kāi)發(fā)等領(lǐng)域。但是由于放射性藥物的低劑量和欠采樣等原因，采集到的投影數(shù)據(jù)經(jīng)常會(huì)被噪聲嚴(yán)重污染，即使最好的圖像重建算法也不能抑制重建圖像中的所有噪聲，并且精確地保護(hù)圖像的邊緣等有用信息，因此尋找一種能夠有效提高重建圖像質(zhì)量的算法是十分必要的。

在經(jīng)典的PET圖像重建算法中，最大似然-期望最大（maximum-likelihood expectation-maximization，MLEM）算法[1]是一種廣泛使用的算法，它的工作原理很簡(jiǎn)單，就是尋找一個(gè)圖像使得測(cè)量數(shù)據(jù)最有可能發(fā)生。和著名的濾波反投影（filtered back-projection，F(xiàn)BP）算法相比，MLEM算法由于考慮了測(cè)量數(shù)據(jù)的統(tǒng)計(jì)特性，其重建圖像的質(zhì)量明顯優(yōu)于FBP算法。但是MLEM算法是一個(gè)“病態(tài)算法”（不適定算法），因?yàn)樗慕庥锌赡懿淮嬖?、不唯一或者不穩(wěn)定[2]?，F(xiàn)在，解決這個(gè)問(wèn)題的有效方法一般是通過(guò)加入一個(gè)正則化項(xiàng)，把不適定的問(wèn)題轉(zhuǎn)化為適定的問(wèn)題，從而達(dá)到提高重建圖像質(zhì)量的目的。這種方法被稱(chēng)為Bayesian算法[3]。Bayesian算法由于具有良好的降噪能力以及解唯一等優(yōu)點(diǎn)，近年來(lái)得到了廣泛的關(guān)注[3-8]。其中，最具代表性的算法是由Green提出的遲一步（one-steplate，OSL）算法[3]，此算法的關(guān)鍵是尋找一個(gè)合適的能量函數(shù)，它直接決定了算法性能的好壞。不幸的是，能量函數(shù)的選擇十分困難。在文獻(xiàn)[4]中，Alenius等人用一個(gè)中值濾波器來(lái)近似代替OSL算法中的能量函數(shù)，提出了著名的中值根先驗(yàn)（median root prior，MRP）算法，此算法雖在一定程度上提高了圖像的質(zhì)量，但是它沒(méi)有與之相對(duì)應(yīng)的目標(biāo)函數(shù)，只是通過(guò)一種經(jīng)驗(yàn)的方式推導(dǎo)出來(lái)，因此不是一種真正意義上的Bayesian算法。2003年，Hsiao等人在MRP算法的基礎(chǔ)上，通過(guò)引入一個(gè)輔助向量，提出了中值先驗(yàn)（median prior,MP）算法[5]，它是一種真正意義上的Bayesian算法。

近年來(lái)，為了進(jìn)一步提高重建圖像的質(zhì)量，基于偏微分方程的圖像重建算法成為新的研究熱點(diǎn)[9-15]。2007年，顏建華和于軍把各向異性擴(kuò)散（anisotropic diffusion，AD）濾波器引入到MP算法中，提出了偏微分方程-中值（partial differential equation median，PDE-median）算法[10]。后來(lái)，桂志國(guó)等人對(duì)PDEmedian算法做了進(jìn)一步的改進(jìn)，如用非局部模糊各向異性擴(kuò)散濾波器代替PDEmedian算法中的AD濾波器等[9]。這些算法都在一定程度上提高了重建圖像的質(zhì)量，但是由于濾波器的缺陷，這些算法對(duì)梯度閾值和擴(kuò)散次數(shù)都十分敏感，實(shí)際運(yùn)用中很難實(shí)現(xiàn)。

為了解決以上算法中參數(shù)設(shè)置困難的問(wèn)題，本文把各向異性中值擴(kuò)散濾波器（anisotropic mediandiffusion，AMD）[16]和MRP算法結(jié)合起來(lái)，提出了一種新的基于Bayesian的PET圖像重建算法，稱(chēng)為各向異性中值擴(kuò)散-中值根先驗(yàn)算法（anisotropic mediandiffusion median root prior，AMDMRP）。新算法能在有效抑制重建圖像中噪聲的同時(shí)精確地保護(hù)圖像的邊緣，并且對(duì)梯度閾值和擴(kuò)散次數(shù)不敏感，易于實(shí)現(xiàn)，實(shí)用性較強(qiáng)。

2　算法基礎(chǔ)

2.1中值根先驗(yàn)算法

1996年，Alenius等人根據(jù)PET圖像的特點(diǎn)，提出了著名的MRP算法。此算法用一個(gè)中值濾波器來(lái)近似代替OSL算法中能量函數(shù)的導(dǎo)數(shù)，是OSL算法的一個(gè)應(yīng)用推廣。MRP算法可用以下公式來(lái)表示[4]：

其中，f(?)表示重建圖像；i表示像素編號(hào)；k表示迭代次數(shù)；表示重建圖像中，編號(hào)為i的像素的鄰域內(nèi)的中值（鄰域的大小一般是一個(gè)3×3的矩形空間）；β是Bayesian參數(shù)，主要用來(lái)設(shè)置先驗(yàn)信息的權(quán)重。

MRP算法的每次迭代過(guò)程都需要計(jì)算MLEM和Bayesian兩個(gè)系數(shù)，MLEM系數(shù)主要負(fù)責(zé)重建圖像，而B(niǎo)ayesian系數(shù)主要負(fù)責(zé)抑制重建圖像中的噪聲，并且保護(hù)圖像的邊緣等細(xì)小結(jié)構(gòu)[4]。如果重建圖像中含有噪聲（局部不單調(diào)），則的值不相等，那么Bayesian系數(shù)的值就不等于0，于是Bayesian系數(shù)就能抑制重建圖像的噪聲。相反，如果重建圖像中不含噪聲（局部光滑），那么由于的值相等，Bayesian系數(shù)的值等于0，從而MRP算法就退化為MLEM算法。MRP算法雖然能在較好地抑制噪聲的同時(shí)保護(hù)重建圖像的邊緣，但是隨著迭代過(guò)程的進(jìn)行，重建圖像中容易產(chǎn)生偽影（“分塊”現(xiàn)象），這些偽影將嚴(yán)重影響臨床診斷的結(jié)果。此外，由于MRP算法沒(méi)有與之相對(duì)應(yīng)的目標(biāo)函數(shù)，只是通過(guò)一種經(jīng)驗(yàn)的方式推導(dǎo)出來(lái)，因而它不是一種真正意義上的Bayesian算法。為了解決這個(gè)問(wèn)題，2003年，Hsiao等人在MRP算法的基礎(chǔ)上，通過(guò)引入一個(gè)輔助向量，提出了一種新的Bayesian算法——中值先驗(yàn)算法（MP）[5]。在PET圖像重建的過(guò)程中，MRP算法的每個(gè)像素值接近于它的局部中值，而MP算法的每個(gè)像素值接近于它的局部中值的中值，因此這兩種算法的重建效果十分相似[5]。

2.2各向異性中值擴(kuò)散濾波器

AD濾波器（又稱(chēng)為P-M模型）是一種非線性、自適應(yīng)的濾波器，目前已經(jīng)廣泛地使用于醫(yī)學(xué)圖像重建等領(lǐng)域，它的基本方程是[17]：

其中，t是一個(gè)時(shí)間參數(shù)；f(x,y,0)表示原始的發(fā)射圖像；div是散度算子；?f(?)是梯度算子；g(?)是擴(kuò)散系數(shù)函數(shù)。Perona和Malik提出了以下兩個(gè)經(jīng)典的擴(kuò)散系數(shù)函數(shù)：

其中，k表示梯度閾值，用來(lái)判斷此區(qū)域是否是一個(gè)局部邊緣。

AD濾波器的基本思想就是根據(jù)圖像內(nèi)容的不同而采取不同的平滑方式，即平滑強(qiáng)度（抑制噪聲的能力）在邊緣區(qū)域減弱，而在非邊緣區(qū)域增強(qiáng)，從而達(dá)到抑制噪聲并且保護(hù)圖像邊緣的目的。因此，P-M模型中梯度閾值的選擇很重要，它直接決定了AD濾波器性能的好壞。但是在一個(gè)被噪聲嚴(yán)重污染的PET圖像中，選取一個(gè)合適的梯度閾值是十分困難的。為了有效地解決這個(gè)問(wèn)題，Ling和Bovid提出了各向異性中值擴(kuò)散濾波器[16]。AMD模型把中值濾波器引入到P-M模型中，每次濾噪的過(guò)程都分為兩步：首先通過(guò)式（2）濾除圖像中具有小梯度值的噪聲；然后再用一個(gè)中值濾波器濾除圖像中具有大梯度值的噪聲。它的擴(kuò)散系數(shù)函數(shù)是：

此外，AMD模型中擴(kuò)散過(guò)程只是濾除具有小梯度值的噪聲，而對(duì)那些具有大梯度值的噪聲，是通過(guò)中值濾波器來(lái)抑制的，因此AMD濾波器對(duì)梯度閾值和擴(kuò)散次數(shù)不敏感。當(dāng)擴(kuò)散次數(shù)不是很小，且梯度閾值不比梯度的標(biāo)準(zhǔn)偏差大很多時(shí)，濾噪效果都很接近[16]。

3　基于AMD模型的Bayesian圖像重建算法

2007年，顏建華和于軍提出的PDEmedian算法在一定程度上提高了重建圖像的質(zhì)量，是一種真正意義上的Bayesian圖像重建算法。在介紹新算法之前，首先簡(jiǎn)單地介紹一下PDEmedian算法的框架。PDEmedian算法是通過(guò)最小化以下兩個(gè)能量函數(shù)得到的[10]：

其中，Ψ(?)是一個(gè)非負(fù)非減函數(shù)；l(f)是對(duì)數(shù)似然函數(shù)；ΦP(f;m)是先驗(yàn)函數(shù)。由于鄰域內(nèi)絕對(duì)值勢(shì)函數(shù)的最優(yōu)化解等于鄰域經(jīng)過(guò)中值濾波的結(jié)果，為了和MRP算法相關(guān)聯(lián)，ΦP(f;m)函數(shù)定義如下：

其中，Nj表示像素 j的鄰域。當(dāng) j′∈Nj時(shí)，wjj′=1，否則，wjj′=0。

接下來(lái)，簡(jiǎn)單分析PDEmedian算法的收斂性。Hsiao等人[5]證明了，如果勢(shì)函數(shù)?是凸的，則先驗(yàn)函數(shù)ΦP(f;m)也是凸的。在式（8）中，?是絕對(duì)值勢(shì)函數(shù)，它是凸的，并且由于l(f)也是凸的，能量函數(shù)E2必是一個(gè)凸函數(shù)。因此，E1和E2的聯(lián)合估計(jì)必定是收斂的，所求的解必定是一個(gè)全局最小解。

PDEmedian算法中梯度閾值K和擴(kuò)散次數(shù)l很重要，它們共同決定了算法性能的好壞。如果K和l的值太小，則不能有效地濾除重建圖像中的噪聲，擴(kuò)散過(guò)程效果不明顯。相反，如果K和l的值太大，則會(huì)使得重建圖像的邊緣模糊，從而丟失大量重要的信息。但是在一個(gè)低信噪比的PET圖像中，噪聲產(chǎn)生的梯度可能超過(guò)邊緣產(chǎn)生的梯度，此時(shí)要選擇一組合適的參數(shù)是十分困難的。因此，對(duì)PDEmedian算法還需要做進(jìn)一步的研究。

本文在PDEmedian算法的基礎(chǔ)上，通過(guò)把AMD模型融合到MRP算法中（由于MP算法與MRP算法的重建結(jié)果十分相似，為了實(shí)驗(yàn)的方便，跟PDEMedian算法一樣，用MRP算法近似代替MP算法），提出了一種新的Bayesian圖像重建算法——AMDMRP算法。新算法的每次迭代更新過(guò)程都可以分為兩步：首先運(yùn)用AMD模型抑制重建圖像中具有小梯度值的噪聲；然后用MRP算法重建圖像，并且進(jìn)一步抑制具有大梯度值的噪聲。具體迭代公式如下：

其中，k表示迭代次數(shù)；j表示像素編號(hào)；f表示重建圖像；Hij是系數(shù)矩陣，表示像素 j發(fā)出的光子，被第i對(duì)探測(cè)器對(duì)捕獲的概率；gi表示第i對(duì)探測(cè)器對(duì)捕獲的光子總數(shù)；w是用來(lái)表示擴(kuò)散率的一個(gè)常量；?fj,j′=fj′-fj；Nj表示 j號(hào)像素的鄰域；Med(?)和β的定義如前所述；擴(kuò)散系數(shù)函數(shù)g(?)的定義如式（5）所示。

考慮到AMD模型對(duì)擴(kuò)散次數(shù)不敏感，為了簡(jiǎn)化AMDMRP算法的參數(shù)設(shè)置，規(guī)定擴(kuò)散次數(shù)和迭代次數(shù)相等，都等于k，那么AMDMRP算法的迭代公式可簡(jiǎn)單地表示為：

AMDMRP算法的具體運(yùn)行步驟如下：

步驟1通過(guò)擴(kuò)散過(guò)程（式（11））抑制噪聲。

步驟2通過(guò)MRP算法重建圖像（式（12））。

步驟3判斷是否滿足迭代終止條件，如果滿足，則轉(zhuǎn)步驟4，否則，轉(zhuǎn)步驟1。

步驟4結(jié)束。

在AMDMRP算法的每次迭代過(guò)程中，擴(kuò)散過(guò)程（式（11））用來(lái)濾除具有小梯度值的噪聲，而具有大梯度值的噪聲卻被接下來(lái)MRP算法（式（12））中的中值濾波器濾除，但是對(duì)那些同樣具有大梯度值的邊緣來(lái)說(shuō)，MRP算法卻不會(huì)影響它們。因此，在AMDMRP算法中，AMD濾波器和MRP算法能夠相互協(xié)同合作，在濾除重建圖像中噪聲的同時(shí)精確地保護(hù)圖像的邊緣。

4　仿真實(shí)驗(yàn)結(jié)果及分析

為了驗(yàn)證AMDMRP算法在PET中的有效性，對(duì)AMDMRP算法和PDEmedian、MRP以及MLEM算法進(jìn)行比較。在模擬仿真實(shí)驗(yàn)中，首先采用的測(cè)試圖像是由計(jì)算機(jī)產(chǎn)生的，大小為128×128像素的Shepp-Logan頭部模型，圖像中每個(gè)像素點(diǎn)的灰度取值范圍都在0～255之間，如圖1所示。假設(shè)投影參數(shù)的大小為128×128，即有128個(gè)投影方向（均勻分布在0～π之間），并且每個(gè)投影方向上有128對(duì)探測(cè)器。利用公式g=Hf產(chǎn)生無(wú)噪聲的觀測(cè)數(shù)據(jù)，用該投影數(shù)據(jù)作為泊松變量的均值生成實(shí)際帶噪聲的投影數(shù)據(jù)。在仿真實(shí)驗(yàn)中，采集到的光子總量約為6×105對(duì)。

Fig.1　Shepp-Logan phantom圖1　Shepp-Logan模型

為了公平，在仿真實(shí)驗(yàn)中所有算法的迭代次數(shù)都設(shè)定為50次，Bayesian參數(shù)β在AMDMRP、PDE-median和MRP算法中都設(shè)置為0.3（這3種算法對(duì)β都不敏感）。在PDEmedian算法中，根據(jù)經(jīng)驗(yàn)（能取得較好的重建結(jié)果）把梯度閾值K和擴(kuò)散次數(shù)l都設(shè)定為40。在AMDMRP算法中，通過(guò)參考圖像的梯度直方圖，把梯度閾值K設(shè)定為10（AMDMRP算法對(duì)梯度閾值不敏感，見(jiàn)圖7）。

4種算法對(duì)頭部模型的重建圖像如圖2所示，圖3是圖2的局部擴(kuò)大圖像。由這些圖可以看出，MLEM算法的重建效果最差，它的重建圖像中含有大量的噪聲；PDEmedian算法的重建效果雖然優(yōu)于MRP和MLEM算法，但是重建圖像中邊緣比較模糊，丟失了大量的細(xì)節(jié)信息。相比之下，AMDMRP算法不僅有效地抑制了重建圖像中的噪聲，還良好地保護(hù)了重建圖像的邊緣，取得了較好的視覺(jué)效果。

Fig.2　Shepp-Logan phantoms reconstructedby 4 different algorithms圖2　4種算法對(duì)Shepp-Logan模型的重建圖像

Fig.3　Zoomed images of Fig.2圖3　圖2的局部擴(kuò)大圖像

為了從客觀上評(píng)價(jià)AMDMRP算法的有效性，使用了NRMSE和SNR兩個(gè)定量評(píng)價(jià)指標(biāo)：

圖4是4種算法重建圖像的NRMSE值隨迭代次數(shù)變化的曲線圖。由圖4可知，新算法的NRMSE值是最小的，說(shuō)明它的重建圖像與原始圖像最接近，重建效果最好。類(lèi)似的結(jié)論也可以通過(guò)分析圖5的SNR變化曲線圖得到。

Fig.4　Plots of NRMSE along with iterations for Shepp-Logan phantoms reconstructed by 4 algorithms圖4　4種算法對(duì)Shepp-Logan模型重建圖像的NRMSE值隨迭代次數(shù)變化的曲線圖

接下來(lái)，通過(guò)與PDEMedian算法的比較來(lái)說(shuō)明AMDMRP算法的穩(wěn)定性（對(duì)梯度閾值和擴(kuò)散次數(shù)不敏感）。在以下實(shí)驗(yàn)中，規(guī)定Bayesian參數(shù)的值β= 0.3。圖6是PDEmedian算法在不同的梯度閾值K和擴(kuò)散次數(shù)l下對(duì)頭部模型的重建圖像，其中迭代次數(shù)k=50。由圖6可知，PDEmedian算法對(duì)梯度閾值和擴(kuò)散次數(shù)十分敏感，它們共同決定了算法重建效果的好壞。當(dāng)K和l的值都比較小時(shí)，算法不能有效地抑制重建圖像中的噪聲，如圖6（a）所示；相反，當(dāng)K 和l的值都比較大時(shí)，雖然能夠抑制重建圖像中的噪聲，但是會(huì)使得圖像的邊緣模糊，從而丟失大量的細(xì)節(jié)信息，如圖6（c）和（d）所示。因此，在PDEmedian算法中選擇一組合適的參數(shù)，使得重建算法既能有效濾噪，又能良好地保護(hù)圖像的邊緣是十分困難的。只有通過(guò)多次反復(fù)的實(shí)驗(yàn)，才能選取一組合適的參數(shù)。在本實(shí)驗(yàn)中，圖6（b）取得了相對(duì)較好的重建效果。圖6中4個(gè)圖的NRMSE和SNR值如表1所示。圖7是AMDMRP算法在不同的梯度閾值K和擴(kuò)散次數(shù)k下對(duì)頭部模型的重建圖像。由圖可知，AMDMRP算法對(duì)梯度閾值不敏感（只要K值不比梯度的標(biāo)準(zhǔn)偏差大很多）。圖7中4個(gè)圖像都很接近，都取得了較好的重建效果。它們的NRMSE和SNR值如表2所示。

Fig.5　Plots of SNR along with iterations for Shepp-Logan phantoms reconstructed by 4 algorithms圖5　4種算法對(duì)Shepp-Logan模型重建圖像的SNR值隨迭代次數(shù)變化的曲線圖

為了進(jìn)一步說(shuō)明算法的有效性，采用胸腔模型為測(cè)試模型。胸腔模型的大小為128×128像素，如圖8所示。投影數(shù)據(jù)的獲取方式和前面的方法類(lèi)似，采集到的光子總數(shù)約為5.4×105對(duì)。

Fig.6　Shepp-Logan phantoms reconstructed by PDEmedian algorithm with different gradient thresholdsKand diffusion numbersl圖6　PDEmedian算法在不同的梯度閾值K和擴(kuò)散次數(shù)l下對(duì)Shepp-Logan模型的重建圖像

Table 1　NRMSE and SNR values for 4 reconstructed images in Fig.6表1圖6中4個(gè)重建圖像的NRMSE和SNR值

Table 2　NRMSE and SNR values for 4 reconstructed images in Fig.7表2圖7中4個(gè)重建圖像的NRMSE和SNR值

Fig.7　Shepp-Logan phantoms reconstructedAMDMRP algorithm with different gradient thresholds Kand diffusion numbersk圖7AMDMRP算法在不同的梯度閾值K和擴(kuò)散次數(shù)k下對(duì)Shepp-Logan模型的重建圖像

同樣，對(duì)AMDMRP算法和PDEmedian、MRP和 MLEM算法進(jìn)行比較。規(guī)定所有算法的迭代次數(shù)都是50，Bayesian參數(shù)β在AMDMRP、PDEmedian和MRP算法中都設(shè)置為0.3。在PDEmedian算法中，根據(jù)經(jīng)驗(yàn)把梯度閾值K和擴(kuò)散次數(shù)l分別設(shè)定為20和40。在AMDMRP算法中，通過(guò)參考圖像的梯度直方圖，把梯度閾值K設(shè)定為20。

Fig.8　Thorax phantom圖8胸腔模型

圖9是以上4種算法對(duì)胸腔模型的重建圖像，它們的局部擴(kuò)大圖像如圖10所示。由這些圖可以看出，AMDMRP算法的重建圖像效果最好，它在濾噪和邊緣保護(hù)兩方面取得了較好的折中。圖11和圖12分別是以上4種算法的NRMSE和SNR值隨迭代次數(shù)變化的曲線圖，它們也客觀地說(shuō)明了AMDMRP算法的有效性。

Fig.9　Thorax phantoms reconstructed by 4 different algorithms圖9　4種算法對(duì)胸腔模型的重建圖像

Fig.10　Zoomed images of Fig.9圖10　圖9的局部擴(kuò)大圖像

圖13是AMDMRP算法在不同梯度閾值和擴(kuò)散次數(shù)下對(duì)胸腔模型的重建圖像。其中，Bayesian參數(shù)的值都是β=0.3。圖中4個(gè)重建圖像十分相似，說(shuō)明了AMDMRP算法對(duì)梯度閾值和擴(kuò)散次數(shù)不敏感。圖13中4個(gè)重建圖像的NRMSE和SNR值如表3所示。

Fig.11　Plots of NRMSE along with iterations for thorax phantoms reconstructed by 4 algorithms圖114種算法對(duì)胸腔模型重建圖像的NRMSE值隨迭代次數(shù)變化的曲線圖

Fig.12　Plots of SNR along with iterations for thorax phantoms reconstructed by 4 algorithms圖12　4種算法對(duì)胸腔模型重建圖像的SNR值隨迭代次數(shù)變化的曲線圖

5　結(jié)束語(yǔ)

為了有效提高PET重建圖像的質(zhì)量，本文通過(guò)把各向異性中值擴(kuò)散濾波器和MRP算法結(jié)合起來(lái)，提出了一種新的Bayesian圖像重建算法——AMDMRP算法。模擬仿真實(shí)驗(yàn)結(jié)果表明，無(wú)論是從主觀的視覺(jué)效果還是客觀的質(zhì)量評(píng)價(jià)參數(shù)來(lái)看，新算法的性能都要好于PDEmedian、MRP和MLEM算法。新算法在抑制噪聲和保護(hù)邊緣兩方面取得了良好的折中，較大程度地提高了重建圖像的質(zhì)量。此外，和PDEmedian算法相比，新算法吸收了各向異性中值擴(kuò)散濾波器的優(yōu)點(diǎn)，在迭代過(guò)程中對(duì)梯度閾值和擴(kuò)散次數(shù)不敏感，易于實(shí)現(xiàn)，實(shí)用性較強(qiáng)。但是還存在一些問(wèn)題，需要予以關(guān)注，主要總結(jié)如下：

Fig.13　Thorax phantoms reconstructed byAMDMRP algorithm with different gradient thresholds Kand diffusion numbersk圖13AMDMRP算法在不同的梯度閾值K和擴(kuò)散次數(shù)k下對(duì)胸腔模型的重建圖像

Table 3　NRMSE and SNR values for 4 reconstructed images in Fig.13表3圖13中4個(gè)重建圖像的NRMSE和SNR值

（1）算法的收斂速度是判斷一個(gè)算法性能好壞的重要指標(biāo)。本文算法沒(méi)有考慮收斂速度的問(wèn)題，下一步可以考慮在收斂速度方面對(duì)它做進(jìn)一步的優(yōu)化，如可以考慮把有序子集的概念用于此算法等。

（2）統(tǒng)計(jì)迭代算法中迭代終止條件的設(shè)置是一個(gè)至今尚未解決的難題，本文是通過(guò)設(shè)置固定迭代次數(shù)來(lái)實(shí)現(xiàn)終止迭代的。下一步思考如何尋找一個(gè)更合適的迭代終止條件來(lái)避免“過(guò)分迭代”的問(wèn)題。

（3）本文算法主要用于二維圖像重建，下一步將考慮三維圖像重建的問(wèn)題。所謂三維圖像重建是指利用二維投影恢復(fù)物體三維信息（形狀等）的數(shù)學(xué)過(guò)程和計(jì)算機(jī)技術(shù)，包括數(shù)據(jù)獲取、預(yù)處理、特征分析和立體匹配等步驟。其中立體匹配是三維圖像重建的核心，它要根據(jù)所提取的特征來(lái)建立圖像對(duì)之間的一種對(duì)應(yīng)關(guān)系，也就是將同一物理空間點(diǎn)在兩幅不同圖像中的成像點(diǎn)一一對(duì)應(yīng)起來(lái)。在進(jìn)行匹配時(shí)要注意場(chǎng)景中一些因素的干擾，比如光照條件、噪聲干擾、物體幾何形狀畸變、表面物理特性以及攝像機(jī)機(jī)特性等諸多變化因素。此外，三維圖像重建的計(jì)算量比二維圖像重建大很多，因此如何提高重建算法的運(yùn)算速度也是三維圖像重建中一個(gè)重要的研究課題。

References:

[1]Shepp LA,Vardl Y.Maximum likelihood reconstruction for emission tomography[J].IEEE Transactions on Medical Imaging,1982,1(2):113-122.

[2]Gaitanis A,Kontaxakis G,Spyrou G,et al.PET image reconstruction:a stopping rule for the MLEM algorithm based on properties of the updating coefficients[J].Computerized Medical Imaging and Graphics,2010,34(2):131-141.

[3]Green P J.Bayesian reconstructions from emission tomography data using a modified EM algorithm[J].IEEE Transactions on Medical Imaging,1990,9(1):84-93.

[4]Alenius S,Ruotsalainen U.Bayesian image reconstruction for emission tomography based on median root prior[J]. European Journal of Nuclear Medicine,1997,24(3):258-265.

[5]Hsiao I T,Rangarajan A,Gindi G.A new convex edgepreserving median prior with applications to tomography[J]. IEEE Transactions on Medical Imaging,2003,22(5):580-585.

[6]Guo Yu,Feng Yuanming,Sun Jian,et al.Automatic lung tumor segmentation on PET/CT images using fuzzy Markov random field model[J].Computational&Mathematical Methods in Medicine,2014(1):171-188.

[7]Norberg J,Roininen L,Vierinen J,et al.Ionospheric tomography in Bayesian framework with Gaussian Markov random field priors[J].Radio Science,2015,50(2):138-152.

[8]Jia Kun,Liang Shunlin,Wei Xiangqin,et al.Automatic land-cover update approach integrating iterative training sample selection and a Markov random field model[J].Remote Sensing Letters,2014,5(2):148-156.

[9]Gui Zhiguo,Liu Yi,He Jiawei.PML algorithm for positron emission tomography combined with nonlocal fuzzy anisotropic diffusion filtering[J].IEEE Transactions on Nuclear Science,2012,59(5):1984-1989.

[10]Yan Jianhua,Yu Jun.Median-prior tomography reconstruction combined with nonlinear anisotropic diffusion filtering[J]. Journal of the Optical Society of America A,2007,24(4): 1026-1033.

[11]He Qian,Huang Lihong.Penalized maximum likelihood algorithm for positron emission tomography by using anisotropic median-diffusion[J].Mathematical Problems in Engineering,2014(1):1-7.

[12]Meng Xianyong,Che Lei,Liu Zhihui,et al.Towards a partial differential equation remote sensing image method based on adaptive degradation diffusion parameter[J].Multimedia Tools&Applications,2015:1-17.doi:10.1007/ s11042-015-2881-1.

[13]Wang Liping,Zhou Shangbo,KarimA.Super-resolution image reconstruction method using homotopy regularization[J]. Multimedia Tools&Applications,2015:1-24.doi:10.1007/ s11042-015-2910-0.

[14]Wang L,Sixou B,Peyrin F.Binary tomography reconstructions with stochastic level-set methods[J].IEEE Signal Processing Letters,2015,22(7):920-924.

[15]Fathi A,Kallivokas L F,Poursartip B.Full-waveform inversion in three-dimensional PML-truncated elastic media[J]. Computer Methods in Applied Mechanics&Engineering, 2015,296:39-72.

[16]Ling J,Bovik A C.Smoothing low-SNR molecular images via anisotropic median-diffusion[J].IEEE Transactions on Medical Imaging,2002,21(4):377-384.

[17]Perona P,Malik J.Scale-space and edge detection using anisotropic diffusion[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1990,12(7):629-639.

HE Qian was born in 1983.He received the M.S.degree from College of Information Engineering,Xiangtan University in 2009.Now he is a Ph.D.candidate at Hunan University,and lecturer at Hunan City University.His research interests include image processing and medical image reconstruction,etc.

何騫（1983—），男，湖南益陽(yáng)人，2009年于湘潭大學(xué)信息工程學(xué)院獲得碩士學(xué)位，現(xiàn)為湖南大學(xué)數(shù)學(xué)與計(jì)量經(jīng)濟(jì)學(xué)院博士研究生，湖南城市學(xué)院講師，主要研究領(lǐng)域?yàn)閳D像處理，醫(yī)學(xué)圖像重建等。

HUANG Lihong was born in 1963.He received the Ph.D.degree from Department of Applied Mathematics,Hunan University in 1996.Now he is a professor and Ph.D.supervisor at Hunan University.His research interests include neural network,differential equation and dynamical system,etc.

黃立宏（1963—），男，湖南湘陰人，1996年于湖南大學(xué)應(yīng)用數(shù)學(xué)系獲得博士學(xué)位，現(xiàn)為長(zhǎng)沙理工大學(xué)黨委委員、副校長(zhǎng)，湖南大學(xué)教授、博士生導(dǎo)師，湖南省數(shù)學(xué)學(xué)會(huì)副理事長(zhǎng)，主要研究領(lǐng)域?yàn)樯窠?jīng)網(wǎng)絡(luò)，微分方程，動(dòng)力系統(tǒng)等。發(fā)表學(xué)術(shù)論文200余篇，主持承擔(dān)973前期研究專(zhuān)項(xiàng)課題1項(xiàng)，國(guó)家自然科學(xué)基金項(xiàng)目6項(xiàng)，其他國(guó)家與省部級(jí)科研和教研項(xiàng)目30余項(xiàng)。

*The National Natural Science Foundation of China under Grant No.11371127(國(guó)家自然科學(xué)基金);the Scientific Research Foundation of Hunan Provincial Education Department under Grant No.15C0253(湖南省教育廳科研項(xiàng)目).

Received 2016-02,Accepted 2016-05.

CNKI網(wǎng)絡(luò)優(yōu)先出版:2016-05-13,http://www.cnki.net/kcms/detail/11.5602.TP.20160513.1434.002.html

文獻(xiàn)標(biāo)志碼：A

中圖分類(lèi)號(hào)：TP391

doi:10.3778/j.issn.1673-9418.1602044

PET Image ReconstructionAlgorithm Combined withAnisotropic Median-Diffusion?

HE Qian1,2+,HUANG Lihong2
1.College of Information Science and Engineering,Hunan City University,Yiyang,Hunan 413002,China
2.College of Mathematics and Econometrics,Hunan University,Changsha 410082,China
+Corresponding author:E-mail:heqian0808@163.com

Abstract:For improving the quality of positron emission tomography(PET)images,this paper proposes a new Bayesian image reconstruction algorithm by combining anisotropic median-diffusion filter with median root prior algorithm.Iterations of the proposed method can be divided into two steps:firstly,suppressing noise with the anisotropic median-diffusion filter;secondly,reconstructing image with median root prior algorithm.Simulation experiment results present that the proposed algorithm can effectively suppress noise and accurately preserve edges information in PET image reconstruction.Furthermore,in comparison to other similar reconstruction algorithms,the proposed method absorbs the advantages of the anisotropic median-diffusion filter and is less sensitive to the selection of the image gradient threshold and diffusion number,thus making the application of PET image reconstruction feasible.

Key words:positron emission tomography(PET);anisotropic median-diffusion;median root prior;image reconstruction;suppressing noise