融合MRI 信息的PET 圖像去噪: 基于圖小波的方法

易利群 盛玉霞 柴 利

正電子發射斷層成像(Positron emission to-mography,PET)是一種強大的功能性成像模式,其通過注入特定的放射性示蹤劑來觀察組織內部的分子水平活動.PET 圖像可以顯示分子水平的示蹤劑分布,從而可以對生理和藥理過程進行定量測量.但是,由于各種物理退化因素和檢測到的光子數量有限,PET 圖像空間分辨率低、噪聲程度高,影響了疾病檢測和診斷的準確性,因此有必要提高PET圖像質量.

對PET 圖像后處理可提升圖像質量,圖像去噪和去模糊是常見的PET 圖像后處理方法.PET圖像的噪聲復雜,沒有精確的模型,通常認為具有泊松和高斯混合噪聲的特性[1],直接應用傳統的去噪方法不能得到好的去噪效果.高斯濾波是一種簡單的去噪算法,常用于PET 圖像去噪,但其在去噪過程中會平滑重要的圖像結構信息.許多學者提出了針對PET 圖像的去噪方法,包括雙邊濾波[2]、各向異性傳播和非局部均值(Non-local means,NLM)濾波[1,3-4]、引導濾波[5]等.文獻[1]提出了一種改進的基于NLM 的動態PET 圖像去噪方法;文獻[4] 提出了一種新的結合多尺度小波變換和NLM 的濾波方法;文獻[5]研究了一種基于磁共振成像(Magnetic resonance imaging,MRI)引導的腦部PET 圖像濾波方法;對PET 圖像去模糊問題,Song 等[6]提出了結合MR (Magnetic resonance)圖像聯合熵先驗的PET 圖像超分辨率重建方法,相對于傳統的全變分和二次罰函數等方法具有更好的去模糊效果;文獻[7]提出了一種基于生成對抗網絡的自監督PET 超分辨率重建方法;稀疏表示字典學習也常用于提升PET 圖像質量[8-9].Li 等[8]通過稀疏表示對低分辨率PET 圖像進行復原;Wang等[9]應用半監督字典學習由低劑量PET 圖像來預測標準劑量圖像;文獻[10]通過多層典型相關分析(Canonical correlation analysis,CCA)方法來估計標準劑量PET 圖像;全變分方法應用于直接重建模型中,以改善PET 圖像重建質量[11-12].Zhang 等[11]提出了一種共邊緣加權全變分正則化的PET-MRI聯合重建方法,該方法可以同時提高兩類圖像的重建質量;文獻[12]提出了一種基于時空全變分正則化的動態PET 圖像重建方法,重建效果優于傳統的最大似然期望最大化(Maximum likelihood expectation maximization,MLEM)方法.以上方法都有效地提高了PET 圖像空間分辨率.

隨著PET/MR 一體化成像設備的應用,高分辨率的MRI 先驗信息可用來進一步提高PET 圖像去噪效果.Yan 等[5]研究了基于MRI 的引導濾波方法,該方法給出結構MR 圖像和功能性PET 圖像之間的局部線性模型,并將部分體積效應(Partial volume correction,PVC)融入到模型中,同時對PET 圖像去噪和部分體積校正.Cheng 等[13]提出了一種將MRI 與迭代高級約束的局部反投影方法(Iterative highly constrained back-projection local region,IHYPR-LR)相結合的去噪方法,該方法將MR 圖像作為IHYPR-LR 方法進行迭代的初始值,在去噪的同時增強了圖像的結構和邊緣.文獻[14]針對腦部成像問題提出了兩種不同的PETMR 核方法來修復PET 圖像.文獻[15]綜述了機器學習在PET 圖像處理中的應用.

近年來,深度學習廣泛應用于醫學圖像處理中[16-18],也為PET 圖像處理提供了新的研究方向[19].文獻[16-17]詳細綜述了深度學習在醫學影像中的最新進展;文獻[19]綜述了基于深度學習的PET 圖像重建方法;文獻[20]提出了基于預訓練VGG (Visual Geometry Group)深度網絡的PET圖像去噪方法,使用人工生成的數據訓練網絡,在實際的PET 圖像上效果一般;文獻[21]通過深度學習方法研究不同級別的低劑量PET 圖像的先驗信息,用來估計高質量的PET 圖像.最近,學者們將生成對抗網絡應用于PET 圖像處理[7,22-23].Song等[7]提出了基于雙重生成對抗網絡的自監督PET超分辨率重建算法,無需配對訓練數據,具有更廣泛的適用性.Gong 等[22]使用混合2D 和3D 網絡學習圖像先驗,復原低劑量PET 圖像.文獻[23]提出了用于PET 圖像去噪的參數傳遞Wasserstein 生成對抗網絡,在保持圖像保真度的同時可以有效地抑制噪聲.Wang 等[24]提出了一種新穎的基于3D自動上下文的局部自適應多模態融合生成對抗網絡模型,用來由低劑量的PET 圖像和MR 圖像合成高質量的PET 圖像,該方法優于深度網絡中傳統的多模態融合方法和PET 估計方法.但醫學圖像存在樣本少且標注困難等問題,深度學習方法在醫學圖像上泛化性能差、網絡結構缺乏通用性、網絡訓練成本大、計算復雜度高.非深度學習方法在樣本少的情況下有其特有的優勢,值得進一步研究.

圖信號處理(Graph signal processing,GSP)是近年興起的高維非規則化數據分析處理新方法,其借助代數圖論和譜圖理論來處理高維加權圖上的信號[25-27].在圖像處理領域,可以通過建立非局部圖像像素間的圖連接關系將圖像轉化為圖信號,圖包含的結構信息可以很好地保留圖像邊緣特性.圖信號處理在圖像濾波、壓縮和復原等應用中取得了較好的效果[28].文獻[29]設計了一種低通圖濾波器,在圖傅里葉域對PET 圖像濾波,去噪效果相比于高斯濾波和NLM 都有所提高,但該方法不能清晰地保留病灶點.Hammond 等[30-31]提出了圖小波變換(Graph wavelet transform,GWT)方法,并給出了GWT 的切比雪夫多項式逼近快速實現算法,文獻[31]詳細介紹了圖小波變換的基礎理論和快速實現方法.相對于圖傅里葉變換(Graph Fourier transform,GFT),圖小波引入多尺度分析,在一些去噪應用中獲得了好的效果[30-32].目前還沒有應用圖小波來對PET 圖像進行去噪的研究.

如上所述,PET 圖像空間分辨率低、噪聲大.已有的PET 去噪方法不能在去噪的同時很好地保留病灶信息.如何既保留PET 圖像的病灶點,又能有效濾除噪聲是PET 圖像去噪的難點問題.針對此問題,本文提出了一種新的基于圖小波的PET圖像去噪方法.圖濾波的去噪性能依賴于圖加權鄰接矩陣的構造,直接用低分辨率PET 圖像構造鄰接矩陣噪聲較大,會影響去噪效果.本文給出了結合MRI 結構信息和PET 圖像的鄰接矩陣構造方法,既能反映MR 圖像的結構信息,又能突出PET圖像的病灶特點,提高圖濾波性能.該方法首先將多幀正弦圖相加后重建,得到PET 合成圖像,然后將PET 合成圖像和MR 圖像融合,由融合圖像構造鄰接矩陣,再對動態PET 圖像進行圖小波變換去噪.仿真數據實驗表明,本文方法與單獨的圖濾波、圖小波變換,以及其他基于模型的PET 圖像去噪方法相比具有更好的去噪性能,與VGG 深度神經網絡等基于學習的方法去噪效果相當,但不需要大量的訓練數據和參數調整.

本文的主要創新性包括: 1)提出了新的基于圖小波的PET 圖像去噪方法;2) 給出了一種融合MRI 結構信息和PET 圖像的構圖方法;3)本文方法能更好地去噪和保持病灶信息,優于圖濾波方法和結合MRI 的傳統PET 圖像去噪方法;4)與深度學習方法相比,本文方法不受標注樣本限制,無需訓練和迭代,計算復雜度低.

1 理論基礎

1.1 圖譜理論

在圖信號處理中,信號x=[x1,x2,···,xN]T通常被轉換為一個在無向加權圖G={V,ε,A}上的信號;其中V={v1,v2,···,vN}是圖上頂點的集合,ε?V ×VT是一組邊,A∈RN×N為加權鄰接矩陣.如果有一條連接了頂點i和頂點j的邊ε=(i,j),則Ai,j ＞0 表示這條邊的權重;否則Ai,j=0.得到鄰接矩陣A后,可以定義圖拉普拉斯矩陣L

其中,D是對角度矩陣,D的對角元素dii是矩陣A第i行元素之和.

圖傅里葉變換將傳統的傅里葉變換推廣到圖譜域,給出了圖信號頂點域與圖譜域之間的傅里葉變換關系.將圖拉普拉斯矩陣分解成特征值和特征向量,令λ?和χ?=[χ?(0),χ?(1),···,χ?(N-1)]T,?=0,1,···,N -1分別為圖拉普拉斯矩陣L對應的特征值和特征向量.特征值表示頻率,特征向量形成圖傅里葉變換基,任何圖信號都可以通過圖拉普拉斯矩陣的特征向量表示.令f=[f(1),f(2),···,f(N)]T為一個圖信號,則它的第i個元素f(i) 是圖G的第i個頂點上的信號,信號f的圖傅里葉變換定義為[25]

其中,f代表頂點域信號,其對應的圖傅里葉逆變換為

1.2 圖小波變換

圖小波是定義在加權圖頂點域數據上的小波變換.圖小波變換需要定義一個圖小波核函數,通過在圖譜域的伸縮操作構造不同的圖譜域小波簇,并在圖譜域調制信號.在文獻[30]中,作者定義圖小波函數算子為Tg=g(L),信號f的圖小波變換在圖譜域的表達式為

其中,g是圖小波核函數,滿足g(0)=0 和limx→∞g(x)=0的條件,它是圖譜域的帶通濾波器,類似于傳統小波變換母函數的傅里葉變換,根據圖傅里葉逆變換(3)得到圖小波逆變換

給定信號f,其圖小波系數和尺度系數定義如下[30]:

1.3 圖小波變換逆變換

圖小波變換通過切比雪夫多項式逼近得到圖小波近似系數,其可以看作是將一個大小為N的輸入信號f映射為大小為N(J+1) 的系數c,即c=W f,其中f∈RN×1為原始信號,J為小波尺度數目.W為矩陣算子,且W存在很多左逆M,使得MW f=f.已知系數c,通過求解方程 (W*W)f=W*c可得到W的偽逆 (W*W)-1W*[30],并將求解出的偽逆作為GWT 的逆[31].為了加快計算速度,采用共軛梯度算法,并應用切比雪夫多項式逼近方法得到重建信號.關于圖小波更詳細的介紹,可參考文獻[30-31].

2 方法

2.1 融合MRI 信息的動態PET 圖像去噪方法

將多個PET 正弦圖短時間幀相加后重建得到合成圖像,可以改善計數統計信息并減少噪聲[33].合成幀圖像盡管丟失了時間信息,但不同時間幀的空間信息可以很好地保留在合成幀中,用于提高短時間幀的重建質量.在本文中,對于有病灶點的PET 圖像,將多幀正弦圖相加得到合成的正弦圖,再通過最大似然期望最大化(Maximum likelihood expectation maximization,MLEM)方法[34]重建得到PET 合成圖像Ic,以保留病灶點信息.為了保留圖像更多的結構信息,本文提出了一種新的結合MRI 信息的融合圖像構造方法,在此基礎上對PET 圖像進行圖小波變換,用來去除PET 圖像中的噪聲.通過硬閾值方法將MR 圖像與PET 合成圖像進行融合,融合方法為

其中,Ic(m,n),IMRI(m,n),IF(m,n) 分別表示PET合成圖像、MR 圖像和融合圖像在點 (m,n) 處的像素值,Threshold為閾值參數.對合成圖像中大于閾值的點,融合圖像取PET 合成圖像的像素值,這部分保留了PET 圖像中的病灶信息;對小于閾值的點,融合圖像取MR 圖像的像素值,以保留結構信息.PET 圖像病灶點的像素值相對于對應位置的MR 圖像像素值偏大,通過設置適當的閾值,可以保留病灶信息.新的融合圖像既保留了PET圖像的病灶信息,也包含有MR 圖像的部分結構信息.

在融合圖像上構造鄰接矩陣A,將融合圖像中每個像素點當作是圖G上的一個頂點,頂點i和頂點j之間的邊緣權重由式(14)的閾值高斯核權重函數定義,即

其中,fin表示PET 圖像Iin的向量形式,切比雪夫多項式階數M=51,j表示圖小波第j個尺度,j=1,2,···,J.對PET 圖像Iin進行第1 次圖小波多尺度分解后,得到表示低頻的圖譜尺度系數1,0和表示高頻的圖小波系數1,j,通過保留圖小波域的低頻段尺度系數1,0,并丟棄高頻段的圖小波系數1,j來達到去噪的效果.同時,為了盡可能保留在高頻段部分丟棄的圖像細節,保留圖小波系數的第1 個尺度系數1,1,進行第2 次圖小波多尺度分解

通過兩次圖小波多尺度分解得到兩組圖小波變換系數,保留兩次分解得到的低頻圖譜尺度系數1,0和2,0.將系數組合=[1,0,2,0,0] 進行圖小波逆變換,通過共軛梯度法求解式(19),得到去噪后的PET 圖像[30]

本文方法的整個實現過程如下:

2.2 參數選擇

本文方法有9 個關鍵參數: 圖像塊的大小n,最近鄰居數目k,平滑參數θ,圖像塊相似度影響參數η,小波尺度數目J,尺度核函數h,小波核函數g,小波尺度參數s和閾值參數Threshold.經過反復實驗,選取相對較優的參數設置,其中尺度核函數h,小波核函數g和小波尺度參數s設置與文獻[30]中相同,其他參數具體設置見表1.

表1 本文方法參數設置Table 1 Parameter setting in this paper

3 實驗與分析

我們在兩個PET 圖像仿真數據集上進行實驗,將本文方法和MRI 引導濾波方法[5]、MRI 引導的IHYPR-LR 方法[13]、直接由PET 合成圖像構圖的圖濾波[29]和圖小波四種方法進行比較.其中圖小波方法是直接由PET 合成圖像構圖,沒有融合MRI信息,相當于閾值為0.其他參數設置與本文融合MRI 信息的方法相同.為公平起見,將對比方法[5,13]的實驗參數都調整到最優.

3.1 評價指標

為了評估和比較各種去噪方法,使用信噪比(Signal-to-noise ratio,SNR)和均方根誤差(Root mean squared error,RMSE) 來評估去噪后的PET 圖像質量.

其中,xi表示原始PET 圖像體素i的強度值,yi表示去噪后的PET 圖像體素i的強度值.

3.2 模擬數據生成

本文使用3D 大腦模型通過計算機模擬Brain-Web 的幻像[35],該模型由脊椎液、白質和灰質組成,并選擇兩個基于18F-FDG 的葡萄糖代謝的組織隔室模型作為動力學模型.本文從正常MR 圖像中選擇了第90 個切片,從異常MR 圖像中選擇了第90個切片,并采用文獻[36]的方法產生區域時間活動曲線(Time activity curve,TAC),在注射示蹤劑時開始動態PET 掃描.動態PET 數據由60 min內的24 個時間幀組成: 4×20 s,4×40 s,4×60 s,4×180 s 和 8×300 s.為了生成帶噪聲的動態PET圖像,本文采用PET 成像領域常用的方法[1],首先通過正向投影生成原始的動態PET 正弦圖,然后將60 min 內產生 7×107個光子總數的泊松過程應用于原始的無噪聲正弦圖中,最后通過60 次MLEM 方法迭代重建得到無病灶動態PET 圖像和單病灶動態PET 圖像.實驗中PET 圖像大小為175×175,MR 圖像大小為 1 81×217,通過Brain-Web 中的配準參數得到已配準的MR 圖像.

3.3 實驗結果與分析

3.3.1 無病灶數據實驗結果與分析

第1 個實驗從光子數為 7×107時生成的無病灶動態PET 大腦圖像中,選擇第6 幀、第12 幀、第18 幀和第24 幀進行去噪效果演示.MR 圖像如圖1(a),從動態PET 圖像中得到的合成圖像如圖1(b)所示,該PET 圖像沒有病灶點,直接使用MR 圖像構圖.觀察圖1(a)可知,MR 圖像含有更多的結構信息,可以幫助我們對動態PET 圖像進行去噪.無病灶情況下PET 圖像去噪實驗結果如圖2 所示,無病灶情況下結合MRI 的PET 圖像去噪方法比較結果見表2,直接由PET 合成圖像構圖的方法與本文方法比較結果見表3.由圖2 可以看出,與其他方法相比,本文方法得到的PET 去噪圖像具有更好的視覺效果.同時在客觀評估指標方面,與結合MRI 的PET 圖像去噪方法、直接由 PET合成圖像構圖的去噪方法相比,本文方法的SNR值最大,RMSE 值最小,說明有更好的去噪效果.

圖1 兩種圖像示例Fig.1 Example of two kinds of images

圖2 無病灶PET 圖像去噪結果Fig.2 Denoising results of normal PET images

表2 無病灶情況下結合MRI 的PET 圖像去噪方法比較Table 2 Comparison of PET image denoising methods incorporated with MRI on the normal dataset

表3 無病灶情況下由PET 合成圖像構圖的方法與本文方法比較Table 3 Comparison of the methods of constructing graph by PET composite image with the proposed method on the normal dataset

3.3.2 單病灶數據實驗結果與分析

第2 個實驗從光子數為 7×107時生成的單病灶動態PET 大腦圖像中,選擇第6 幀、第12 幀、第18 幀和第24 幀進行去噪效果演示.BrainWeb中相應的MR 圖像如圖3(a)所示,動態PET 合成圖像如圖3(b)所示,由MR 圖像和PET 合成圖像得到的融合圖像如圖3(c)所示.觀察圖3(c)可以看出,融合圖像既保留了MRI 圖像的結構信息,也保留了合成圖像中的病灶點信息.單病灶情況下5種方法去噪后的PET 圖像如圖4 所示,由圖4 可以看出,本文方法在去噪的同時還保留了更多的結構信息和病灶信息,有更好的視覺效果.各種方法的SNR 和RMSE 值如表4 和表5 所示.可以看出,本文方法的SNR 值最大,RMSE 值最小.為進一步評估本文方法的去噪性能,我們分別畫出了噪聲圖像和5 種方法去噪后的PET 圖像殘差圖,如圖5所示.觀察圖5 可知,本文方法得到的PET 圖像與原圖更相近,并且病灶點和邊緣細節保留得更完整.

圖3 單病灶圖像Fig.3 Single-lesion images

圖4 單病灶PET 圖像去噪結果Fig.4 Denoising results of single-lesion PET images

圖5 單病灶PET 圖像去噪殘差圖Fig.5 Denoising residual map of single-lesion PET images

表4 單病灶情況下結合MRI 的PET 圖像去噪方法比較Table 4 Comparison of PET image denoising methods incorporated with MRI on the single-lesion dataset

表5 單病灶情況下由PET 合成圖像構圖的方法與本文方法比較Table 5 Comparison of the methods of constructing graph by PET composite image with the proposed method on the single-lesion dataset

為了驗證不同噪聲級別下本文方法的有效性,在單病灶情況下,本文設置了光子數分別為7×108、7×109兩種情況來模擬生成不同噪聲級別的PET 圖像,各種方法的SNR 和RMSE 值如表6 所示.可以看出,本文方法在不同噪聲程度下,SNR值最大,RMSE 值最小,有更好的去噪效果.

表6 單病灶PET 圖像在不同光子數時各種去噪方法比較Table 6 Comparison of different denoising methods for single-lesion PET images with different photon numbers

為了定量比較不同方法在病灶點區域的重建效果,用病灶點區域對比度恢復系數(Contrast recovery coefficient,CRC)與白質區域(背景區域)標準差(Standard deviation,STD)作為衡量指標.CRC 越接近1,說明病灶點區域的重建去噪效果越好;STD 值越小,說明背景區域的噪聲越少.CRC與STD 的計算式為[20]

本文對比高斯(Gaussian)濾波、MRI 引導濾波(GuideMRI)[5]、局部反投影方法(HYPRMRI)[13]、圖小波(GWT)和本文方法(GWTMRI)對病灶點的恢復效果.指標參數的設置與文獻[20]相同.用以上5 種方法對MLEM 重建圖像去噪,MLEM 共迭代120 次,每迭代24 次計算一個CRC-STD 值,改變迭代次數得到的CRC-STD 曲線如圖6 所示.可以看出,本文方法相對于高斯濾波、GuideMRI、HYPRMRI、GWT 有更好的CRC-STD 折中值.在同樣水平的STD 下,本文方法有更大的CRC 值.因沒有算法實現代碼,不能直接與基于深度學習的方法對比,我們間接比較各種方法的CRC-STD 圖.從文獻[20]中的深度學習方法與高斯濾波的CRCSTD 曲線對比圖可以看出,文獻[20]中的方法和本文方法相對于高斯濾波性能都有很大程度的提升.值得指出的是,該仿真示例僅表明了本文所提方法具有與基于深度學習方法相當的去噪性能.但本文方法的優勢在于可解釋性好,不需要大規模數據訓練.

圖6 不同方法病灶點CRC 與背景區域STD 曲線圖Fig.6 CRC-STD curves of the different denoising methods for single-lesion PET images

4 結束語

提出了一種融合MRI 信息的圖小波變換方法用于動態PET 圖像去噪.本文將MR 圖像與PET合成圖像通過硬閾值方法進行融合,在融合圖像上構造圖拉普拉斯矩陣,最后通過圖小波變換對動態PET 圖像去噪.在兩種不同的BrainWeb 仿真數據上進行實驗,實驗結果表明,與結合MRI 的PET圖像去噪方法以及基于圖濾波的去噪方法相比,本文方法在圖像視覺效果和評估指標方面的表現都是最優的.這說明本文提出的融合MRI 信息的圖小波方法可以更大程度地降低噪聲,同時保留更多的PET 圖像病灶信息,具有更好的去噪性能.本文使用的圖小波函數是文獻中給定的常用函數,根據PET 圖像的特點,能否設計性能更好的專門針對醫學圖像的圖小波函數,值得進一步研究.