甄琰明，張金銘，曲桂紅，袁 波，趙書俊

[1.鄭州大學(xué)物理學(xué)院(微電子學(xué)院)，河南 鄭州 450001；2.北京大基康明醫(yī)療設(shè)備有限公司，北京 100176]

PET是核醫(yī)學(xué)領(lǐng)域中比較先進的成像技術(shù)，其圖像重建算法可分為解析法和迭代法。針對不同算法對PET圖像質(zhì)量的影響已有較多研究[1]。濾波反投影(filter back projection， FBP)算法[2]重建圖像偽影較多，相比有序子集期望最大化(ordered subsets expectation maximization，OSEM)算法圖像質(zhì)量較差[3-4]。臨床PET圖像重建多采用OSEM算法[5]。由于放射性藥物低劑量和欠采樣等原因，采集到的投影數(shù)據(jù)(正弦圖)經(jīng)常被噪聲嚴重污染而影響重建圖像質(zhì)量[6]。為提高重建圖像質(zhì)量，可在圖像重建后、正弦圖中或基于迭代統(tǒng)計的重建過程中進行噪聲平滑。在正弦圖中進行濾波處理計算效率高，無需過多修改已有機器配置，并能在重建圖像中產(chǎn)生均勻的各向同性采樣分辨率。MOKRI等[7-8]提出3D均值中值濾波器能有效去除噪聲，并保持圖像的邊緣，但并未討論濾波參數(shù)選擇及不同濾波參數(shù)對重建圖像的影響。本研究根據(jù)投影數(shù)據(jù)梯度分布直方圖確定濾波參數(shù)的選取范圍，從中等間隔地均勻選取一組濾波參數(shù)，對投影數(shù)據(jù)進行3D均值中值濾波，利用開源斷層圖像重建(software for tomographic image reconstruction， STIR)軟件[9-10]中的3D OSEM算法重建濾波后投影數(shù)據(jù)，以視覺和定量方法評估濾波參數(shù)對3D OSEM重建圖像的影響。

1 原理與方法

1.1 仿真數(shù)據(jù) 采用ECAT 962 PET掃描儀，利用分析模擬軟件(ASIM)[11]模擬維度為288(radial bins)×144(angular views)×239(slices)的3D PET模體數(shù)據(jù)[12]，并根據(jù)米歇爾圖(Michelogram)對其進行組合排列[13]，對包括600萬真實符合事件、600萬隨機符合事件和210萬散射符合事件在內(nèi)的數(shù)據(jù)進行過衰減校正和軸向歸一化校正。

1.2 3D OSEM算法 OSEM算法按投影角度將投影數(shù)據(jù)分成有限個有序子集，每個子集應(yīng)用最大似然期望最大化(maximum likelihood expectation maximization， MLEM)[14]算法對圖像更新一次,為一次子迭代,所有子集全部遍歷為一次完整的迭代。相對MLEM算法,OSEM算法在一次迭代過程中重建圖像更新n次，圖像收斂速度加快[15-16]。OSEM的迭代公式如下：

(1)

式中yi表示第i條射線的投影數(shù)據(jù)，Hij表示第j號體素發(fā)出的光子被第i對探測器對探測到的概率，f表示待重建圖像，k為該算法的迭代次數(shù)，Ts為第s個子集，s=1,2,...,n。

1.3 3D均值中值濾波算法[7-8]均值中值濾波器是基于最小化角度模糊偽影、平滑平坦區(qū)域及保持圖像邊緣3個目標(biāo)設(shè)計的。定義3D正弦圖為I(l,θ,z)，l是徑向投影，θ是角度，z是切片。對正弦圖的每個點定義一個梯度參數(shù)h(h∈[0,1])，根據(jù)高斯濾波后的正弦圖的局部梯度信息得到該參數(shù)，計算公式如下：

(2)

F=(f2l+f2θ+f2z)1/2

(3)

式中fl、fθ和fz分別是高斯濾波后的正弦圖在l、θ和z方向上的一階導(dǎo)數(shù)。圖像存在邊緣時，一定有較大的梯度值；相反，圖像比較平滑部分灰度值變化較小，則相應(yīng)梯度也較小。h表示歸一化后的梯度值，用于控制每個正弦圖點上的平滑和增強水平，在平坦區(qū)域處應(yīng)用較多平滑，在邊緣區(qū)域處應(yīng)用較少平滑。根據(jù)h值選取1個閾值Tm，將正弦圖劃分為2個區(qū)域進行分區(qū)域濾波[8]，見圖1。

圖1 區(qū)域1和區(qū)域2的劃分

對區(qū)域1中的數(shù)據(jù)用2種不同尺寸模板計算均值，濾波后區(qū)域1中的所有數(shù)據(jù)值將替換為：

(4)

MeanA(l,θ,z)和MeanB(l,θ,z)分別表示整個3D數(shù)據(jù)用模板A和模板B完成均值濾波后的值。模板A比模板B具有更高平滑度。公式(4)表示平滑隨h值增加而逐漸減少。h=0時，只采用模板A來施加大量平滑；h∈[0,Tm]之間，模板A的平滑程度通過結(jié)合模板B的貢獻而減慢；h最接近Tm點處時，平滑程度則由模板B控制，即對較小h值施加大量平滑，對較大h值施加少量平滑，以去除噪聲保持邊緣。

圖2 梯度分布直方圖

對區(qū)域2中數(shù)據(jù)采用角度分別為0°、15°、45°、75°、90°、105°、135°、165°的8個方向模板計算均值中值。方向模板表示正弦圖中點的局部性，通過包括比角度分量更多的徑向分量來減少重建圖像中的角度模糊偽影。濾波后區(qū)域2中的所有數(shù)據(jù)值將替換為：

IRegion2(l,θ,z)=(1-h)×Meandir(l,θ,z)+h×Mediandir(l,θ,z)

(5)

均值Meandir(l,θ,z)和中值Mediandir(l,θ,z)的計算同理，中值計算公式如下:

(6)

ωi=ω1×ω2

(7)

mi是數(shù)組M的第i個方向模板Md的中值，為每個mi分配相應(yīng)的權(quán)重ωi。如公式(8)計算權(quán)重ω1，σstd是標(biāo)準(zhǔn)差。如公式(9)計算權(quán)重ω2，dr是每組模板Md的平均值與中值的絕對差值，而σr是總計數(shù)值的平均值與中值的絕對差值。

(8)

(9)

對區(qū)域2中的濾波采用基于h值的加權(quán)策略，而不選擇具有最佳權(quán)重的單個模板來估計均值Meandir(l,θ,z)和中值Mediandir(l,θ,z)，目的在于確保正弦圖強度值隨強度降低而逐漸平穩(wěn)變化，避免正弦圖中強度突然變化導(dǎo)致重建圖像產(chǎn)生切向條紋，特別是在高強度區(qū)域和低強度區(qū)域之間的邊界處。

2 結(jié)果

2.1 實驗結(jié)果 用MATLAB畫出正弦圖梯度分布直方圖(圖2)，據(jù)以確定梯度控制參數(shù)K值的選擇范圍，從其中等間隔地均勻取K值分別為100、150、200、300、400、500、600。由公式(2)可知，參數(shù)h∈(0,1)，在0到1之間等間隔選取閾值Tm。分別用不同參數(shù)對數(shù)據(jù)進行均值中值濾波處理，得到數(shù)組數(shù)據(jù)。利用STIR中的3D OSEM算法對濾波前后3D正弦圖數(shù)據(jù)進行重建，設(shè)置子集個數(shù)為9，迭代次數(shù)為12，得到數(shù)組實驗結(jié)果。重建圖像噪聲及邊緣保持效果對閾值Tm并不敏感，而對濾波參數(shù)K十分敏感；且實驗結(jié)果符合MOKRI等[8]提出的高計數(shù)正弦圖需要較低Tm值的結(jié)論，因此選取閾值較小的數(shù)組圖像(圖3～5)進行視覺評估,著重分析濾波參數(shù)K對重建圖像質(zhì)量和效果的影響。

參數(shù)K值過小時，圖像邊緣保持良好，但內(nèi)部細節(jié)模糊；若參數(shù)K值過大，則邊緣缺失，引起圖像過光滑。不論如何選取K值，對比未濾波重建圖像，濾波后重建圖像中的噪聲光圈消失，表明濾波后噪聲減少，達到了去噪的目的。參數(shù)K值處于最大與最小中間的某個值時，既能有效抑制噪聲，又可保持良好的圖像邊緣，且對比度有所提高。

2.2 定量分析 利用信噪比(signal noise ratio, SNR)、標(biāo)準(zhǔn)差(standard deviation, STD)和Brenner梯度函數(shù)定量評價仿真重建圖像。圖6中ROI局部SNR、局部STD計算結(jié)果見表1、2。K值≥300時，ROI的SNR一直小于未濾波數(shù)據(jù)重建圖像SNR，若要保持良好圖像邊緣，則此區(qū)間范圍的值不可取。濾波后重建圖像ROI的STD恒小于未濾波數(shù)據(jù)重建圖像ROI的STD，表明濾波后重建圖像噪聲減小，達到去噪目的。

局部SNR定義為ROI的局部均值與局部STD之比，如公式(10)所示。

(10)

Brenner梯度函數(shù)是最常用的計算圖像清晰度的梯度函數(shù)[17]。正焦的清晰圖像比模糊的離焦圖像邊緣更加銳利清晰，邊緣像素灰度值變化大，因而梯度值更大，故可利用梯度函數(shù)獲取圖像灰度信息評判圖像清晰度?；贐renner梯度函數(shù)的圖像清晰度定義如公式(11)所示。

圖3 采用3D OSEM重建濾波前后數(shù)據(jù)，原始數(shù)據(jù)重建圖像(A)及固定Tm為0.1及濾波參數(shù)K分別為100、200、600的重建圖像(B～D)圖4 采用3D OSEM重建濾波前后數(shù)據(jù)，原始數(shù)據(jù)重建圖像(A)及固定Tm為0.2及濾波參數(shù)K分別為100、200、600的重建圖像(B～D) 圖5 采用3D OSEM重建濾波前后數(shù)據(jù)，原始數(shù)據(jù)重建圖像(A)及固定Tm為0.3及濾波參數(shù)K分別為100、200、600的重建圖像(B～D)

表1 不同濾波參數(shù)邊緣局部SNR和STD

注：未濾波數(shù)據(jù)重建圖像的邊緣局部SNR=3.63，局部STD=0.17

表2 不同濾波參數(shù)中央局部SNR和STD

注：未濾波數(shù)據(jù)重建圖像的中央局部SNR=3.36，局部STD=0.81

(11)

其中I(x,y)表示圖像I對應(yīng)像素點(x,y)的灰度值，D(I)為圖像清晰度計算結(jié)果。

計算結(jié)果見表3。K值≥300時，濾波后數(shù)據(jù)重建圖像的清晰度一直小于未濾波數(shù)據(jù)重建圖像，提示處于此范圍的值不可取。在考慮背景噪聲減少和邊緣信號保留及圖像清晰度的情況下，綜合圖像視覺效果和定量分析結(jié)果，在梯度分布直方圖區(qū)間100～300之間折中選取一個K值。本研究中，K值為150時實現(xiàn)了去除噪聲保持邊緣目的，可提供視覺清晰、高分辨率的圖像。

不同數(shù)據(jù)的梯度分布直方圖的梯度值范圍不同。為使實驗結(jié)果更具有參考價值，計算梯度分布直方圖中K<150時的梯度分布占比為96.19%，作為參考標(biāo)準(zhǔn)。K<150時的梯度分布占比=K<150時的梯度計數(shù)值/K<600時的梯度計數(shù)值。對于不同數(shù)據(jù)，若求得梯度分布占比在96.19%左右，可自濾波參數(shù)范圍中確定一個適當(dāng)濾波參數(shù)進行折中，以獲得視覺效果清晰的高質(zhì)量圖像。

圖6 ROI A.局部SNR； B.局部STD

表3 不同濾波參數(shù)圖像清晰度D(I)

注：未濾波數(shù)據(jù)重建圖像的清晰度D(I)=1.36

3 討論

PET可在病變區(qū)域內(nèi)物理結(jié)構(gòu)變化尚不明顯的早期通過發(fā)現(xiàn)新陳代謝異常而加以確定，從而實現(xiàn)疾病早發(fā)現(xiàn)、早診斷及早治療，提高治愈率。為臨床診斷提供視覺清晰的高分辨率圖像一直是重要研究方向之一。

本研究采用的均值中值濾波器的基本思想是根據(jù)不同h內(nèi)容采取不同平滑方式，平滑強度(抑制噪聲能力)在邊緣區(qū)域減弱，在非邊緣區(qū)域增強，從而抑制噪聲，保護圖像邊緣。參數(shù)K值直接影響h中梯度信息的分布，梯度控制參數(shù)K直接決定均值中值濾波器的性能。

本研究結(jié)果表明，選取K值直接影響重建圖像的噪聲及邊緣保持效果。圖像中存在邊緣部分時一定有較大梯度值，相反，圖像中有比較平滑部分時，灰度值變化較小，相應(yīng)梯度也較小。如所選K值過小，會將包含邊緣信息的部分大梯度值劃分到區(qū)域2，對區(qū)域2數(shù)據(jù)同時進行均值中值濾波可更好地保留邊緣信號，但不能抑制圖像中的所有噪聲，尤其是具有大梯度值的噪聲；相反，選取K值過大易將包含邊緣信息的部分大梯度值劃分到區(qū)域1，而區(qū)域1數(shù)據(jù)只接受均值濾波，雖可顯著消除噪聲，但會導(dǎo)致圖像邊緣缺失，使圖像過度平滑。因此，K值最好是大于噪聲產(chǎn)生的梯度并小于邊緣產(chǎn)生的梯度。定量分析結(jié)果顯示，濾波后數(shù)據(jù)重建圖像的中央局部STD恒小于濾波前數(shù)據(jù)重建圖像，可知不論如何選取K值，均可達到去噪目的。綜合主觀視覺效果和定量分析結(jié)果，可在區(qū)間100～300內(nèi)確定一個最優(yōu)K值。本研究在此區(qū)間內(nèi)折中選取K值，從視覺效果和定量評估中獲得視覺清晰、高分辨率的圖像，并計算該值在梯度分布直方圖中的梯度分布占比，將其作為一個參考標(biāo)準(zhǔn)。針對不同數(shù)據(jù)可依據(jù)此標(biāo)準(zhǔn)從濾波參數(shù)范圍中確定適當(dāng)?shù)臑V波參數(shù)，達到折中效果，獲得視覺清晰的高質(zhì)量圖像。

本研究所涉方法及數(shù)據(jù)僅針對本機型，且僅初步分析不同濾波參數(shù)對3D OSEM重建圖像的影響，并確定濾波參數(shù)的選取范圍；今后將進一步尋求最優(yōu)化算法來選取最優(yōu)值。