上官宏，劉 祎，張 權(quán)，桂志國＋

（1.中北大學(xué) 電子測試技術(shù)國防重點(diǎn)實(shí)驗(yàn)室，山西 太原030051；2.中北大學(xué) 儀器科學(xué)與動態(tài)測試教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原030051）

0 引 言

正電子發(fā)射斷層成像（positron emission tomography，PET），融合了核物理學(xué)、計(jì)算機(jī)圖形學(xué)、醫(yī)學(xué)影像學(xué)和分子生物學(xué)等前沿學(xué)科，能夠代表現(xiàn)今核醫(yī)學(xué)技術(shù)的先進(jìn)水平，在腦功能成像和癌癥的早期診斷中提供重要的檢測信息［1，2］。Bayesian理論指出，MAP算法在ML類算法的基礎(chǔ)上添加一項(xiàng)正則化約束即先驗(yàn)分布信息，能夠更加準(zhǔn)確的描述投影數(shù)據(jù)的物理模型。但是先驗(yàn)信息選取不當(dāng)時(shí)，則可能影響重建圖像的平滑、邊緣及噪聲特性。所以如何準(zhǔn)確地描述先驗(yàn)信息，對PET 重建具有非常重要的意義［3－5］。1982年，Shepp和Vardi首先提出了MLEM 算法，該算法能夠考慮系統(tǒng)的物理效應(yīng)和噪聲的統(tǒng)計(jì)特性，但速度較慢容易產(chǎn)生棋盤效應(yīng)。2007年，Yan提出了迭代間各向異性擴(kuò)散濾波的正則化MLEM（MLEM－PDE），該算法在去除噪聲的同時(shí)，能很好的保持圖像的邊緣信息。但是MLEM－PDE所采用的P－M（Perona－Malik）擴(kuò)散模型僅考慮圖像的梯度信息，無法有效保持低梯度值處的細(xì)節(jié)信息，且對孤立噪聲的去除效果不佳［5－9］。本文提出了一種結(jié)合圖像灰度方差的正逆各項(xiàng)異性擴(kuò)散模型，該模型采用圖像的梯度算子和局部灰度方差共同表征圖像的特征，同時(shí)考慮前向擴(kuò)散和后向擴(kuò)散，對圖像的不同區(qū)域進(jìn)行不同的處理。此外，在模糊理論的基礎(chǔ)上改進(jìn)了MLEM－PDE 算法，形成了基于灰度方差和正逆擴(kuò)散結(jié)合的優(yōu)質(zhì)PET 重建算法（GVFAB－PML）。實(shí)驗(yàn)結(jié)果表明，新算法克服了傳統(tǒng)模型造成的過度平滑，重建出更加優(yōu)質(zhì)的圖像。

1 基本理論

1.1 PET統(tǒng)計(jì)模型

1982年，Shepp和Vardi首先提出了觀測數(shù)據(jù)的泊松模型，在PET 重建中，采用的離散模型形式為

式中：y＝［y1，y1，…，yN］T表示觀測投影數(shù)據(jù)；f＝［f1，f2，…，fM］T表示圖像向量，H 為系統(tǒng)矩陣，ε為非泊松噪聲?；谧畲笏迫还烙?jì)理論，MLEM 算法表示為

1.2 P－M 擴(kuò)散模型

各向異性擴(kuò)散濾波器源于經(jīng)典的熱傳導(dǎo)方程，Perona和Malik提出的各向異性擴(kuò)散模型如下式所示

這里f是梯度算子，g（ f ）為控制擴(kuò)散強(qiáng)度的擴(kuò)散系數(shù)函數(shù)。1990年，Perona和Malik提出

2002年，Gilboa等提出正逆各向異性擴(kuò)散模型，把圖像中梯度較大的地方定義為邊緣；梯度較低的部分定義為背景和噪聲，對圖像的不同區(qū)域進(jìn)行不同的處理。其表達(dá)式如下

式中：Kf——正向擴(kuò)散閾值與K 作用相同大致相同，Kb——逆向擴(kuò)散中心點(diǎn)，與ω共同控制后向擴(kuò)散的作用范圍，α調(diào)節(jié)正向和逆向擴(kuò)散的比重。

1.3 非局部MRF二次先驗(yàn)?zāi)Ｐ?/h3>

2005年，Buades等最先提出了NL－Means算法，基本思想是：當(dāng)前像素值由圖像中所有與它結(jié)構(gòu)相似的像素加權(quán)平均得到。2009 年，Chen等［10］基于非局部均值思想設(shè)計(jì)了一種非局部MRF二次先驗(yàn)?zāi)Ｐ?，如下?/p>

式中：UNL——先驗(yàn)?zāi)芰浚豮j＿NL——非局部像素k和像素j之間的權(quán)值，Z（k）——配分函數(shù)。

2 重建算法

2.1 結(jié)合圖像灰度方差的正逆各項(xiàng)異性擴(kuò)散模型

針對正逆各項(xiàng)異性擴(kuò)散模型僅考慮圖像梯度信息的不足，本文將2010年Chao［11，12］提出的圖像灰度方差引入擴(kuò)散模型，用圖像灰度方差和梯度共同表征圖像的特征。認(rèn)為邊緣鄰域內(nèi)的梯度值和灰度方差都比較大，而背景區(qū)和含有豐富細(xì)節(jié)處的梯度值比較小，但后者處灰度方差要明顯大于前者。為了防止細(xì)節(jié)信息丟失，提出的改進(jìn)擴(kuò)散模型如下所示

記f 為圖像的灰度值，在3×3鄰域中計(jì)算其灰度方差

式中：珚f——3×3鄰域中的平均灰度。因?yàn)榛叶确讲畹姆秶@著大于梯度的范圍，將灰度方差按比例縮小，使灰度方差和梯度相互兼容。對于一個(gè)8位灰度圖像，標(biāo)準(zhǔn)化方差定義為

2.2 基于灰度方差和正逆擴(kuò)散結(jié)合的優(yōu)質(zhì)PET重建算法

傳統(tǒng)最大似然算法考慮了系統(tǒng)的物理模型和觀測數(shù)據(jù)的統(tǒng)計(jì)特性，重建結(jié)果較FBP 算法更加優(yōu)越。但隨迭代次數(shù)增加，MLEM 棋盤效應(yīng)嚴(yán)重，重建結(jié)果不穩(wěn)定導(dǎo)致圖像退化。針對此問題，本文使用先驗(yàn)懲罰項(xiàng)對MLEM 算法進(jìn)行約束，提出了基于灰度方差和正逆擴(kuò)散結(jié)合的優(yōu)質(zhì)PET重建算法（GVFAB－PML）。每次迭代重建之前，先對當(dāng)前重建圖像進(jìn)行正逆擴(kuò)散濾波，然后將Chen 等［12］提出的MRF先驗(yàn)?zāi)Ｐ鸵氲阶畲笏迫凰惴ㄖ校捎脩土P最大似然算法重建。該算法主要分為兩大步驟：

（1）基于改進(jìn)正逆擴(kuò)散模型的各向異性擴(kuò)散

（2）PML重建算法

若將MRF非局部二次先驗(yàn)?zāi)Ｐ妥鳛閼土P項(xiàng)引入到最大似然算法中，則目標(biāo)函數(shù)變?yōu)?/p>

繼而，得出其迭代更新方程

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

3.1 實(shí) 驗(yàn)

本文采用Matlab2008進(jìn)行算法仿真，利用兩組體模數(shù)據(jù)進(jìn)行重建：體模1為shepp－logan頭部剖面模型；體模2為Thorax胸腔截面模型，體模大小均為128×128，投影數(shù)據(jù)Sinogram 光子總量近似為1.01×106，投影參數(shù)為：128×128，即128個(gè)投影方向（均勻分布在0°～180°內(nèi)）并且每個(gè)方向上有128個(gè)探測器對。系統(tǒng)矩陣A 經(jīng)過歸一化處理，其規(guī)模為（128×128）×（128×128）。在仿真實(shí)驗(yàn)中，算法的參數(shù)都是根據(jù)經(jīng)驗(yàn)和重建結(jié)果手動設(shè)定的最優(yōu)值。本文所有方法迭代次數(shù)為100，全局控制參數(shù)設(shè)為0.9105，鄰域窗口大小為3×3。體模數(shù)據(jù)和觀測數(shù)據(jù)如圖1所示。

本文設(shè)計(jì)了一系列仿真實(shí)驗(yàn)，來驗(yàn)證算法有效性，本文算法GVFAB－PML、PML－PDE 與NL－PML，以及傳統(tǒng)的MAP、MLEM 等的重建結(jié)果如圖2、圖4所示，各迭代算法重建結(jié)果的局部放大圖如圖3所示。比較圖2（a）～圖2（g）和圖4（a）～圖4（g）可以看出，GVFAB－PML、PMLPDE、NL－PML重建效果較傳統(tǒng)算法優(yōu)越性顯著，可見各項(xiàng)異性擴(kuò)散思想的引入極大的改善了重建圖像質(zhì)量。比較圖2（d）～圖2（f）和圖4（d）～圖4（f）可知GVFAB－PML、PML－PDE、NL－PML噪聲平滑效果都比較好，其中新算法GVFAB－PML去噪效果最好，重建圖像與真實(shí)圖像很接近。在邊緣保持方面，比較圖3（d）～圖3（f）可以看出，NL－PDE算法中間的細(xì)小物體部分有大量偽影產(chǎn)生邊緣不清晰。PML－PDE算法破壞了物體的邊緣結(jié)構(gòu)，特別是中間的細(xì)小部分產(chǎn)生了明顯的形變。GVFAB－PML算法相對而言能較好的保持邊緣，基本保留原有特征。綜上，新算法GVFAB－PML在噪聲去除和邊緣保持方面表現(xiàn)良好，具有優(yōu)越性。

圖1 體模數(shù)據(jù)和觀測數(shù)據(jù)

圖2 體模1各重建算法結(jié)果

圖3 體模1各重建算法結(jié)果局部放大

圖4 體模2各重建算法結(jié)果

3.2 重建結(jié)果量化分析

本文采用均方誤差（root mean squared error，RMSE）和信噪比（signal to noise ratio，SNR）來量化分析各重建算法的性能

從圖5（a）可知，本文算法GVFAB－PML具有最小的均方誤差，這說明它的重建結(jié)果與原始圖像最接近。且隨著迭代次數(shù)的增加，曲線趨于平穩(wěn)。同理從圖5（b）可以看出，GVFAB－PML算法信噪比最大，在多次迭代后仍然具有較好的噪聲特性，逐漸收斂趨于平穩(wěn)，算法優(yōu)越性得到直觀體現(xiàn)。為了更加清楚的描述各質(zhì)量參數(shù)，表1給出了本文所用重建算法的量化數(shù)據(jù)。

分析表1可知，由于各項(xiàng)異性擴(kuò)散思想的引入，GVFAB－PML、PML－PDE算法較NL－PDE，以及傳統(tǒng)的MAP、MLEM 算法在噪聲平滑方面有大幅度的改善。與其它幾種算法相比，NL－PML 算法的去噪表現(xiàn)非常優(yōu)異，量化的反映了非局部均值濾波和各向異性擴(kuò)散算法在降噪方面的優(yōu)越性。同理分析均方誤差，MAP算法均方誤差值很大，說明邊緣細(xì)節(jié)信息丟失，偽影嚴(yán)重。GVFAB－PML、PML－PDE算法較NL－PML，以及傳統(tǒng)的MAP、MLEM 的重建質(zhì)量得到改善，細(xì)節(jié)更加清晰。從表中數(shù)據(jù)可以看出，本文算法GVFAB－PML 在細(xì)節(jié)保持方面的改善最為明顯。綜上分析，GVFAB－PML算法在圖像質(zhì)量改善和噪聲平滑方面都有改進(jìn)，算法可行性得到驗(yàn)證。

圖5 各重建算法評估參數(shù)變化曲線

表1 圖1中所示重建圖像的評估對比

4 結(jié)束語

本文提出了一種結(jié)合圖像灰度方差的正逆各項(xiàng)異性擴(kuò)散模型，該模型采用圖像的梯度算子和局部灰度方差共同表征圖像的特征，同時(shí)考慮正向擴(kuò)散和逆向擴(kuò)散，對圖像的不同區(qū)域進(jìn)行不同的處理。文中GVFAB－PML 算法，結(jié)合了模糊理論和各向異性擴(kuò)散的優(yōu)點(diǎn)，克服了傳統(tǒng)模型造成的過度平滑，保留了細(xì)節(jié)信息。

設(shè)計(jì)了一系列仿真實(shí)驗(yàn)，來驗(yàn)證所提算法的優(yōu)越性。給出了整體和局部對比圖，對比分析了各算法的噪聲平滑效果及細(xì)節(jié)特性。實(shí)驗(yàn)結(jié)果表明，無論從視覺效果還是客觀準(zhǔn)則來看，GVFAB－PML 算法在圖像質(zhì)量改善和噪聲平滑方面都有改進(jìn)，算法可行有效。

［1］Rahmim A，Zaidi H.PET versus SPECT：Strengths，limitations and challenges ［J］.Nuclear Medicine Communications，2008，29 （3）：193－207.

［2］Rahmim A，Dinelle K，Cheng J C，et al.Accurate event－driven motion compensation in high－resolution pet incorporating scattered and random events［J］.IEEE Transactions on Medical Imaging，2008，27 （8）：1018－1033.

［3］ZHOU Jian.Image reconstruction for positron emission tomography ［D］.Nanjing：Southeast University，2006 （in Chi－nese）.［周健.正電子發(fā)射斷層的圖像重建方法研究 ［D］.南京：東南大學(xué)，2006.］

［4］HUANG Jing，MA Jianhua，LU Lijun，et al.Generalized Gibbs priors in positron emission tomography ［J］.Dianzi Xuebao（Acta Electronica Sinica），2010，38 （4）：899－903 （in Chinese）.［黃靜，馬建華，路利軍，等.基于廣義Gibbs先驗(yàn)的優(yōu)質(zhì)PET 成像 ［J］.電子學(xué)報(bào)，2010，38 （4）：899－903.］

［5］YAN Jianhua.Investigation of positron emission tomography image reconstruction ［D］.Wuhan：Huazhong University of Science ＆Technology，2007 （in Chinese）.［顏建華.正電子發(fā)射斷層圖像重建算法研究 ［D］.武漢：華中科技大學(xué)，2007.］

［6］HE Jiawei.The research of positron emission tomography statistical iterative algorithm and acceleration methods ［D］.Taiyuan：North University of China，2011 （in Chinese）.［何佳偉.正電子發(fā)射斷層成像的統(tǒng)計(jì)迭代方法及加速方法研究［D］.太原：中北大學(xué)，2011.］

［7］Gui Z，Liu Y，He J.PML algorithm for positron emission tomography combined with nonsocial fuzzy anisotropic diffusion filtering ［J］.IEEE Transactions on Nuclear Science，2012，59 （5）：1984－1989.

［8］Gui Z，He J.Regularized maximum likelihood algorithm for PET image reconstruction using a detail and edges preserving anisotropic diffusion ［J］.Optik－International Journal for Light and Electron Optics，2012，123 （6）：507－510.

［9］Zhang Q，Liu Y，Shu H，et al.Application of regularized maximum likelihood algorithm in PET image reconstruction combined with nonlocal fuzzy anisotropic diffusion ［J］.Optik－International Journal for Light and Electron Optics，2013，124（20）：4561－4565.

［10］CHEN Yang.High quality PET image reconstruction based on different prior acquisitions［D］.Guangzhou：Southern Medical University，2007 （in Chinese）.［陳陽.基于不同先驗(yàn)獲取的PET 圖像優(yōu)質(zhì)重建新方法研究 ［D］.廣州：南方醫(yī)科大學(xué)大學(xué)，2007.］

［11］Chao S M，Tsai D M.An improved anisotropic diffusion model for detail－and edge－preserving smoothing ［J］.Pattern Recognition Letters，2010，31 （13）：2012－2023.

［12］Chao S M，Tsai D M.Anisotropic diffusion with generalized diffusion coefficient function for defect detection in low－contrast surface images ［J］.Pattern Recognition，2010，43 （5）：1917－1931.