高小帆，張 權(quán)，劉 祎，張 芳，孫未雅，桂志國(guó)，2

（1.中北大學(xué) 電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，山西 太原030051；

2.中北大學(xué) 儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原030051）

0 引 言

隨著計(jì)算機(jī)斷層掃描（CT）技術(shù)在醫(yī)學(xué)早期預(yù)測(cè)與防治中的廣泛應(yīng)用，以及公眾自我健康意識(shí)的提高，在不影響診斷質(zhì)量的前提下，人們開(kāi)始追求盡可能減少放射劑量.但由于放射劑量的降低會(huì)使重建圖像噪聲增大，質(zhì)量退化，繼而嚴(yán)重影響醫(yī)學(xué)診斷與治療，因此在降低輻射劑量時(shí)，重建出高質(zhì)量圖像的研究備受關(guān)注［1］.Rust［2］等使用非線性高斯濾波器鏈對(duì)重建的圖像進(jìn)行濾波，得到了不錯(cuò)的降噪效果，很好地保持了邊緣和細(xì)節(jié).Lui D［3］提出了一種新穎的噪聲補(bǔ)償CT重建方法，提高了重建圖像的信噪比.王麗艷［4］等利用待重建圖像稀疏性的先驗(yàn)信息作為正則項(xiàng)，泊松噪聲的負(fù)對(duì)數(shù)似然函數(shù)作為保真項(xiàng)來(lái)設(shè)計(jì)優(yōu)化目標(biāo)函數(shù)，從而達(dá)到去噪且保持細(xì)節(jié)的目的.Chen Yang［5］等通過(guò)使用非局部自適應(yīng)加權(quán)先驗(yàn)統(tǒng)計(jì)方法對(duì)圖像進(jìn)行重建，改善了低劑量CT圖像的質(zhì)量.

低劑量CT投影數(shù)據(jù)噪聲模型的數(shù)據(jù)特點(diǎn)已被廣泛研究，其中Wang J［6］等人通過(guò)對(duì)多個(gè)體模反復(fù)進(jìn)行實(shí)驗(yàn)和分析，得出低劑量CT投影數(shù)據(jù)經(jīng)過(guò)對(duì)數(shù)變換后，其均值和方差呈非線性遞增關(guān)系，近似服從非平穩(wěn)高斯分布的結(jié)論.由于最大似然期望最大算法（MLEM）在重建過(guò)程中考慮了觀測(cè)數(shù)據(jù)的統(tǒng)計(jì)特性，且其具有非負(fù)性、全局收斂性和計(jì)數(shù)保持的特點(diǎn)，能較好地改善重建圖像的效果，該算法被廣泛運(yùn)用于對(duì)投影數(shù)據(jù)統(tǒng)計(jì)特性進(jìn)行分析.在實(shí)際運(yùn)用中，迭代次數(shù)達(dá)到一定數(shù)量后，隨著迭代次數(shù)的增加重建圖像的質(zhì)量會(huì)出現(xiàn)棋盤(pán)效應(yīng)，產(chǎn)生失真及非收斂的迭代過(guò)程.基于Bayesian理論的最大后驗(yàn)（Maximum A Posterior，MAP）方法，有效地解決了此問(wèn)題［7］，該算法考慮了低劑量CT投影數(shù)據(jù)的統(tǒng)計(jì)特性，通過(guò)對(duì)先驗(yàn)分布加入先驗(yàn)信息，使重建過(guò)程中進(jìn)行多次迭代后仍然很好地抑制了噪聲，克服了MLEM重建算法收斂慢的缺點(diǎn).MAP重建的思想主要是在傳統(tǒng)的MLEM算法基礎(chǔ)上加入先驗(yàn)約束，從而達(dá)到抑制噪聲、平滑圖像、增強(qiáng)圖像邊緣的目的［8］.但是基于傳統(tǒng)的貝葉斯方法提供的先驗(yàn)信息是有限的，往往會(huì)使低劑量重建圖像出現(xiàn)階梯狀偽影和過(guò)平滑現(xiàn)象［9］.本文將從此問(wèn)題出發(fā)對(duì)低劑量CT重建進(jìn)行進(jìn)一步的研究.

基于偏微分方程（Partial Differential Equation，PDE）各向異性擴(kuò)散的降噪算法［10-12］，可以滿足圖像不同強(qiáng)度去噪的需求，即是一種自適應(yīng)的去噪技術(shù)，近年來(lái)被眾多學(xué)者廣泛地運(yùn)用.傳統(tǒng)的二階偏微分降噪會(huì)出現(xiàn)“階梯”效應(yīng)，而四階偏微分可以有效抑制此缺點(diǎn)，既可以克服“階梯”效應(yīng)，又可以根據(jù)梯度和切線方向的不同擴(kuò)散程度，有效地保持邊緣.本文受文獻(xiàn)［13］的啟發(fā)把該文獻(xiàn)中提出的四階偏微分算法用于本文，并且在此基礎(chǔ)上與絕對(duì)差值排序檢測(cè)法進(jìn)行了結(jié)合，對(duì)低劑量CT圖像進(jìn)行降噪處理；同時(shí)在每次迭代過(guò)程中，把降噪后的重建圖像引入MP重建算法的目標(biāo)函數(shù)，提出一種基于改進(jìn)四階各向異性擴(kuò)散的中值先驗(yàn)重建算法，得到了很好的重建效果.實(shí)驗(yàn)結(jié)果表明，本文算法既改善了圖像的質(zhì)量，又提高了圖像的信噪比.

1 MP重建算法

Alenius等在1997年提出了MRP（Median Root Prior）重建方法.MRP法使圖像的像素接近其領(lǐng)域的中值［14］.該算法能有效地保持邊緣，但不是實(shí)際意義的MAP算法，只是一個(gè)經(jīng)驗(yàn)公式.2003年，Hsiao等在MRP的基礎(chǔ)上，構(gòu)造了輔助向量，從而提出基于中值先驗(yàn)（Median Prior，MP）的重建方法［15］，此算法是真正意義上的MAP法.

Hsicao定義的中值先驗(yàn)?zāi)繕?biāo)函數(shù)公式為

式中：Φ（y|f）為對(duì)數(shù)型似然函數(shù)；R（f，m）為一種新型先驗(yàn)分布目標(biāo)函數(shù)；y表示觀測(cè)數(shù)據(jù)向量；f表示圖像向量；m是輔助向量，與f具有相同的維數(shù).

先驗(yàn)分布的目標(biāo)函數(shù)如下式所示：

式中：φ為先驗(yàn)勢(shì)函數(shù)；Nj是像素j的鄰域；ωjk為權(quán)值因子，表示fj和輔助向量鄰域像素mf′的相互作用的強(qiáng)弱，當(dāng)j∈Nj時(shí)，ωjk=1，否則ωjk=0.MP使用的是絕對(duì)值勢(shì)函數(shù)：R（f，m）=輔助向量鄰域像素mf′選取為mj′=median｛fj，w｝，w表示所選中值算子的窗口大小.

2 改進(jìn)的重建算法

2.1 基于絕對(duì)差值排序檢測(cè)的四階偏微分方程降噪

基于偏微分方程（Partial Differential Equation，PDE）的各向異性擴(kuò)散降噪法，屬于自適應(yīng)降噪技術(shù)，該技術(shù)是在圖像的平滑區(qū)域增強(qiáng)平滑強(qiáng)度，而在邊緣區(qū)域則自適應(yīng)地削弱平滑強(qiáng)度，從而在降噪的同時(shí)避免破壞圖像的邊緣.由于二階PDE降噪會(huì)出現(xiàn)“階梯”效應(yīng)，文獻(xiàn)［13］中提出了一種四階PDE降噪算法，具體的分析以及完整的公式詳見(jiàn)文獻(xiàn)，本文給出主要公式為

式中：fηη和fξξ分別為圖像梯度方向和切線方向的二階導(dǎo)數(shù)，表達(dá)式為為擴(kuò)散系數(shù)函數(shù).由擴(kuò)散系數(shù)函數(shù)表達(dá)式可知：在圖像的邊緣處，梯度值較大，函數(shù)的值較??；在圖像的平坦區(qū)域，梯度值較小，函數(shù)的值較大.由此可知，擴(kuò)散系數(shù)函數(shù)利用圖像梯度信息作為邊緣檢測(cè)算子來(lái)控制圖像的平滑程度.

式（3）的離散化形式為

四階偏微分降噪算法雖然可以克服二階偏微分降噪出現(xiàn)的階梯現(xiàn)象，但其去噪能力相對(duì)較弱，如果為了提升效果而增加迭代次數(shù)，又會(huì)對(duì)圖像的邊緣造成破壞.Roman Garnett等人在文獻(xiàn)［16］中提出了一個(gè)局部圖像統(tǒng)計(jì)——絕對(duì)差值排序 檢 測(cè) 法（rank-ordered absolute differences，ROAD），該方法可以表示出當(dāng)前像素值與其周圍鄰域像素值強(qiáng)度的差異.其具體原理在文獻(xiàn)［16］中有詳細(xì)的闡述，這里不再贅述，本文受該方法啟發(fā)，用于四階各向異性擴(kuò)散的擴(kuò)散函數(shù)中，其主要公式為

式中：dx，y表示像素x與y強(qiáng)度差的絕對(duì)值，并將d，中的值按升序進(jìn)行排列，定義

式中：ri（x）表示第i小的dx，y.

ROAD提供了一個(gè)計(jì)算當(dāng)前像素與它周圍鄰域像素之間相似度的方法.考慮到圖像的內(nèi)部區(qū)域和邊緣的連續(xù)性，故在像素周圍的八鄰域中，至少有四個(gè)鄰域值與當(dāng)前像素的強(qiáng)度相似，說(shuō)明ROAD值較小.而噪聲會(huì)使當(dāng)前像素與大多數(shù)鄰域像素的強(qiáng)度相差較大，說(shuō)明ROAD值較大.因此，通過(guò)ROAD可以用于區(qū)分邊緣與噪聲，故把ROAD用于四階各向異性擴(kuò)散的擴(kuò)散函數(shù)中是合適的.

2.2 基于改進(jìn)的四階各向異性擴(kuò)散的MP重建算法

由于MAP方法引入先驗(yàn)信息，從而改善了低劑量CT重建圖像的質(zhì)量與提高圖像的信噪比，在一定程度上保持了圖像的細(xì)節(jié)信息［17］.基于傳統(tǒng)的貝葉斯法只能提供有限的先驗(yàn)信息，故使圖像過(guò)平滑且出現(xiàn)偽影.而改進(jìn)的四階偏微分對(duì)噪聲的敏感性較高，可以根據(jù)梯度和切線方向不同的擴(kuò)散程度，從而有效地保持邊緣.故把改進(jìn)的四階偏微分處理和MAP算法相結(jié)合是合適的.由于基于傳統(tǒng)的貝葉斯方法只能提供有限的局部先驗(yàn)信息，故把改進(jìn)的四階偏微分處理的結(jié)果運(yùn)用到MP的目標(biāo)函數(shù)中作為本文的正則項(xiàng)可以得到令人滿意的效果.

本文算法的目標(biāo)函數(shù)為

式中：Φ（y|f）為對(duì)數(shù)型的似然函數(shù)；U（f）為一種新的先驗(yàn)分布的正則模型.這里U（f）取與四階各向異性擴(kuò)散相關(guān)的先驗(yàn).

式中：F（uj）為對(duì)圖像進(jìn)行改進(jìn)的四階偏微分處理；這里φ（Δ）=Δ2／2.

經(jīng)過(guò)計(jì)算，最后的表達(dá)式為

2.3 基于改進(jìn)四階各向異性擴(kuò)散的中值先驗(yàn)重建算法

基于改進(jìn)四階偏微分先驗(yàn)的MP重建算法增加了基于改進(jìn)四階偏微分的正則項(xiàng)，可以有效地改善重建圖像的質(zhì)量，但是依然會(huì)有些許塊狀的偽影，故本文在每次迭代中，對(duì)圖像進(jìn)行基于改進(jìn)四階偏微分先驗(yàn)的MP重建算法處理后，再對(duì)重建的圖像進(jìn)行改進(jìn)的四階偏微分降噪.實(shí)驗(yàn)結(jié)果表明，本文的算法在降噪的同時(shí)可以很好地保持圖像的細(xì)節(jié)和邊緣信息，提高圖像的信噪比.

本文的具體重建算法如下：

1）基于改進(jìn)的四階偏微分先驗(yàn)的MP重建算法，見(jiàn)式（11）.

2）在對(duì)圖像進(jìn)行基于四階偏微分先驗(yàn)的MP重建算法處理后，對(duì)重建的圖像進(jìn)行四階偏微分降噪

其中，擴(kuò)散函數(shù)為

3）重復(fù)以上過(guò)程一定次數(shù)后，得到最終的重建圖像.

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

3.1 重建圖像比較

本文所有算法的實(shí)驗(yàn)仿真環(huán)境均為：計(jì)算機(jī)配置為Windows 7旗艦版32位SP1（Direct X 11）的操作系統(tǒng)，處理器是英特爾Celeron（賽揚(yáng)）E3300＠2.50 GHz雙核，內(nèi)存為2 GB.編程工具使用MATLAB7.6.0（R2008a）.本文首先選取大小為128 mm×128 mm的Sheep-Logan頭部剖面圖模型作為實(shí)驗(yàn)對(duì)象，如圖1（a）所示.本文選取大小為16 384×16 384的系統(tǒng)矩陣，所有實(shí)驗(yàn)均采用平行投影的方式，在180個(gè)角度中取128個(gè)投影方向，每個(gè)方向分配128對(duì)探測(cè)器.本文在理想投影數(shù)據(jù)中加入如下式關(guān)系的高斯噪聲后對(duì)低劑量的投影數(shù)據(jù)進(jìn)行仿真.

式中：i=1，2，…，N為探測(cè)器信道，N是信道總數(shù)；λi為第i個(gè)探測(cè)器獲得的投影數(shù)據(jù)的平均值；σ2i為第i個(gè)探測(cè)器獲得投影數(shù)據(jù)的方差；ki為第i個(gè)探測(cè)器的參數(shù)；T為系統(tǒng)參數(shù)；對(duì)于給定的CT采集系統(tǒng)，ki與T是給定的.本文參數(shù)選取為：ki=200，T=12 000.為了驗(yàn)證本文算法的有效性，將本文算法與傳統(tǒng)的MLEM，MRP（Median Root Prior），對(duì)圖像的每次迭代中進(jìn)行MLEM重建后直接進(jìn)行四階各向異性擴(kuò)散的基于四階的算法，以及在每次重建迭代中使用基于方差的各向異性擴(kuò)散降噪的算法進(jìn)行了比較.

各種算法中涉及到的各種參數(shù)以及迭代次數(shù)，均為反復(fù)實(shí)驗(yàn)后得到的最優(yōu)值，圖1中的每個(gè)重建圖像都是達(dá)到最優(yōu)時(shí)得到的圖像.由圖可知，傳統(tǒng)的MLEM重建圖像的質(zhì)量最差；MRP算法的結(jié)果圖可以達(dá)到降噪的效果，但是有明顯的塊狀偽影；基于四階和基于方差重建算法對(duì)重建圖像的噪聲進(jìn)行了一定程度的抑制且可獲得比較清晰的圖像，但是圖像中存在比較明顯的塊狀陰影；本文算法在消除噪聲時(shí)保持了圖像的邊緣和細(xì)節(jié)信息，使圖像達(dá)到比較優(yōu)質(zhì)的效果，從視覺(jué)上分析，重建效果和其他幾種算法相比達(dá)到最優(yōu)，初步表明本文算法是有效的.

圖2為胸腔模型的實(shí)驗(yàn)結(jié)果.由圖可知，傳統(tǒng)的MLEM重建圖像的質(zhì)量最差；MRP算法的結(jié)果圖可以達(dá)到降噪的效果，但是有明顯的塊狀偽影；基于四階和基于方差重建算法對(duì)重建圖像的噪聲進(jìn)行了一定程度的抑制且可獲得比較清晰的圖像，但是圖像中存在一些比較明顯的塊狀陰影；本文算法在解決低劑量CT圖像噪聲問(wèn)題的同時(shí)，對(duì)圖像的邊緣和細(xì)節(jié)信息也有較好的保持，且提高了圖像的信噪比，從視覺(jué)上分析，重建效果優(yōu)于其他幾種算法，進(jìn)一步說(shuō)明該算法的有效性.

圖3為骨骼組織模型圖像的實(shí)驗(yàn)結(jié)果.所得結(jié)論與圖1實(shí)驗(yàn)結(jié)論相同，再一次說(shuō)明了該算法的有效性.

圖1 Sheep-Logan頭部剖面模型圖像處理結(jié)果對(duì)比 Fig.1 Comparison of processing results on head profile model of Sheep-Logan

圖2 胸腔模型圖像處理結(jié)果對(duì)比 Fig.2 Comparison of processing results on thoracic model

圖3 骨骼組織模型圖像處理結(jié)果對(duì)比 Fig.3 Comparison of processing results on skeletal tissues model

3.2 重建精度比較

由上述分析可知，本文提出的算法對(duì)低劑量CT的重建圖像有很好的處理能力，既可以保持圖像的邊緣又達(dá)到降噪的目的.本文采用以下指標(biāo)對(duì)重建圖像的質(zhì)量進(jìn)行定量的描述.

1）歸一化均方誤差（Root Mean Squared Error，NMSE）

2）歸一化均方距離（Normalized Mean Square Distance，NMSD）

3）均方絕對(duì)誤差（Mean Absolute Error，MAE）

4）信噪比（Signal to Noise Ratio，SNR）

式中：J表示圖像像素點(diǎn)的總和；Fi和fi分別表示重建圖與原始圖的第個(gè)像素的灰度值；Mi和mi分別表示重建圖與原始圖的均值.這些指標(biāo)從不同的方面評(píng)價(jià)重建圖像與原始圖像的接近程度以及重建圖像的質(zhì)量，表1～表3為本文算法與其他算法的客觀評(píng)價(jià)結(jié)果.

表1 Sheep-Logan頭部剖面圖模型各種算法的客觀評(píng)價(jià) Tab.1 Objective evaluation of different algorithms of the head profile model of Sheep-Logan

表2 胸腔模型各種算法的客觀評(píng)價(jià) Tab.2 Objective Evaluation of different algorithms of thoracic model

表3 骨骼組織模型各種算法的客觀評(píng)價(jià) Tab.3 Objective Evaluation of different algorithms of skeletal tissues

從表1～表3可得出本文算法的信噪比最大，且其他指標(biāo)值均為最小.該結(jié)論說(shuō)明本文算法的重建圖像和原始圖像最為接近.因此在定量評(píng)價(jià)方面，同樣可以表明本算法在低劑量CT重建中是可行的.

圖4，圖5分別給出了本文所用的Sheep-Logan模型和胸腔模型的原始圖與各種算法的重建結(jié)果圖的側(cè)面輪廓線的比較結(jié)果.從圖中可以看出本文方法的重建圖像與原始圖像的吻合度是最高的，最接近于理想圖像，具有最小的噪聲波動(dòng)，故本文算法可以有效地解決低劑量CT重建圖像的噪聲問(wèn)題.

圖4 頭部剖面圖模型各種算法第65行側(cè)面輪廓線的對(duì)比 Fig.4 Comparison of various algorithms of the head profile model on the side of the contour line 65th

圖5 胸腔模型各種算法第65行側(cè)面輪廓線的對(duì)比 Fig.5 Comparison of various algorithms of thoracic model on the side of the contour line 65th

4 結(jié) 論

本文提出了一種基于改進(jìn)的四階各向異性擴(kuò)散的中值先驗(yàn)重建算法.該算法先在中值先驗(yàn)MP算法的基礎(chǔ)上，對(duì)目標(biāo)函數(shù)進(jìn)行了修訂，加入改進(jìn)的四階各向異性擴(kuò)散正則項(xiàng)，該正則項(xiàng)可以對(duì)重建圖像進(jìn)行降噪的同時(shí)保持圖像的邊緣和細(xì)節(jié)，從而形成基于改進(jìn)四階各向異性的MP重建算法；該算法可以很好地改善圖像的質(zhì)量，但是依然會(huì)存在一些塊狀偽影；四階各向異性擴(kuò)散對(duì)噪聲的敏感性比較高，可以根據(jù)梯度和切線方向不同的擴(kuò)散程度，進(jìn)而有效地保持圖像的邊緣，結(jié)合這兩種方法的優(yōu)勢(shì)，本文提出了基于改進(jìn)的四階各向異性擴(kuò)散的中值先驗(yàn)重建算法，即在每次迭代中，對(duì)基于改進(jìn)的四階各向異性擴(kuò)散的MP重建算法處理后的圖像，再進(jìn)行改進(jìn)四階各向異性擴(kuò)散降噪處理，從而進(jìn)一步提高圖像的抗噪聲性能.實(shí)驗(yàn)結(jié)果表明，該算法無(wú)論在主觀效果還是客觀效果上，均說(shuō)明該算法是可行的.

［1］Linton O W，Mettler F A.National conference on dose reduction in CT with an emphasis on pediatric patients［J］.American Journal of Roentgenology，2003，181（2）：321-329.

［2］Rust G F，Aurich V，Reiser M.Noise dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20 m As with nonlinear Gaussian filter chains［C］∥Medical Imaging 2002 Confer-ence，New York：IEEE，2002：186-197.

［3］Lui D，Cameron A，Modhafar A，et al.Low-dose computed tomography via spatially adaptive Monte-Carlo reconstruction［J］.Computerized Medical Imaging and Graphics，2013，37（7-8）：438-449.

［4］王麗艷，韋志輝.低劑量CT的線性Bregman迭代重建算法［J］.電子與信息學(xué)報(bào)，2013，35（10）：2418-2424.Wang Liyan，Wei Zhihui.Linearized bregman iterations for low-dose CT reconstruction［J］.Journal of Electronics and Information，2013，35（10）：2418-2424.（in Chinese）

［5］Chen Y，Gao D Z，Nie C，et al.Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting local nonprior［J］.Computerized Medical Imaging and Graphics，2009，33（7）：495-500.

［6］Li T，Li X，Wang J，et al.Nonlinear sinogram smoothing for low-dose X-ray CT［J］.IEEE Transactions on Nuclear Science，2004，51（5）：2505-2513.

［7］Li S T.Markov random field modeling in image analysis［M］.Berlin：Springer Verlag，2001.

［8］何玲君，潘晉孝，孔慧華.自適應(yīng)正則MAP的CT圖像重建方法研究［J］.計(jì)算機(jī)工程與應(yīng)用，2011，47（28）：198-200.He Lingjun，Pan Jinxiao，Kong Huihua.Adaptive regularized MAP of CT image reconstruction method［J］.Computer Engineering and App Hcations，2011，47（28）：198-200.（in Chinese）

［9］李曉紅，張權(quán)，劉祎，等.基于小波收縮和正逆擴(kuò)散結(jié)合的Pattern Recognition Letters質(zhì)中值先驗(yàn)圖像重建算法［J］.計(jì)算機(jī)應(yīng)用，2012，32（12）：3357-3360.Li Xiaohong，Zhang Quan，Liu Yi，et al.High quality median prior image reconstruction algorithm based on wavelet shrinkage and forward-and-backward diffusion［J］.Journal of Computer Applications，2012，32（12）：3357-3360.（in Chinese）

［10］Cho S I，Kang S J，Kim H S，et al.Dictionary-based anisotropic diffusion for noise reduction［J］.Pattern Recognition Letters，2014，46：36-45.

［11］Veerakumar T，Esakkirajan S，Vennila I.Edge preserving adaptive anisotropic diffusion filter approach for the suppression of impulse noise in images［J］.AEU-International Journal of Electronics and Communications，2014，68（5）：442-452.

［12］Surya Prasath V B，Vorotnikov D.Weighted and well-balanced anisotropic diffusion scheme for image denoising and restoration［J］.Nonlinear Anal.Real World Appl，2014，17：33-46.

［13］Mohammad Reza Hajiaboli.An anisotropic fourth-order diffusion filter for image noise removal［J］.International Journal of Computer Vision，2011，92（2）：177-191.

［14］Karali E，Koutsouris D.Towards novel regularization approaches to PET image reconstruction［J］.Journal of Biosciences and Medicines，2013，1（2）：6-9.

［15］Hsiao I T，Rangarajan A，Gindi G.A new convex edge-preserving median prior with applications to tomography［J］.IEEE Transactions on Medical Imaging，2003，22（5）：580-585.

［16］Roman Garnett，Timothy Huegerich，Charles Chui.A universal noise removal algorithm with an impulse detector［J］.IEEE Transactions on image processing，2005，14（11）：1747-1754.

［17］Zhan Jie，Chen Wufan.Bayesian reconstruction algorithm for PET using new Markov quadratic hybrid multi-order priors［C］∥The 1st International Conference on Bioinformatics and Biomedical Engineering.New York：IEEE，2007：334-337.