本發(fā)明屬于pet成像技術領域，具體涉及一種基于塊字典學習和稀疏表達的pet圖像重建方法。

背景技術：

正電子發(fā)射斷層掃描(positronemissiontomography，pet)是放射型斷層成像的一種，依靠進入病人體內(nèi)的放射性示蹤劑進行斷層成像。pet圖像重建主要分為解析方法和迭代方法(均為統(tǒng)計方法)，其中解析方法發(fā)展較早，比較著名的是濾波反投影方法(filteredback-projection，fbp)。fbp方法基于雷登變換(radontransform)，因為沒有考慮到儀器在測量時的噪聲，其重建準確率不高。隨著計算機性能的提高，迭代方法成為發(fā)展熱點，但是迭代方法如極大似然-期望值最大化(maximumlikelihood-expectationmaximization，ml-em)方法在迭代次數(shù)高的時候存在病態(tài)解等問題。加入懲罰項的極大后驗估計方法(mostaposterior，map)考慮了圖像的先驗知識，一定程度上克服了ml-em方法存在的問題，于是懲罰項的設計成了學者關注的焦點?；谙∈璞磉_和字典學習的方法因為能夠借助解剖圖像如ct和mri的先驗信息，因此可以設計成懲罰項。另一方面，塊字典學習相比于普通的字典學習方法，能夠獲得更小的表示誤差和更快的匯聚速度。

技術實現(xiàn)要素：

鑒于上述，本發(fā)明提出了一種基于塊字典學習和稀疏表達的pet圖像重建方法，即利用塊字典學習方法學習字典和稀疏表達作為懲罰項(先驗函數(shù))，利用基于泊松分布假設的似然函數(shù)，形成一個最大后驗估計問題。

一種基于塊字典學習和稀疏表達的pet圖像重建方法，包括如下步驟：

(1)利用探測器對注入有放射性示蹤劑的生物組織進行探測，根據(jù)各個位置探測器對所采集到的符合計數(shù)向量，構建得到符合計數(shù)矩陣y；

(2)根據(jù)pet成像原理建立pet的測量方程如下，通過對所述測量方程引入泊松噪聲約束，得到pet基于泊松分布的似然函數(shù)l(x)；

y＝gx+r+v

其中：g為系統(tǒng)矩陣，x為pet濃度分布矩陣，r和v分別為關于反射符合事件和散射符合事件的測量噪聲矩陣；

(3)利用解剖圖像通過訓練學習得到塊字典d，并根據(jù)其稀疏表達形成pet基于塊字典的懲罰函數(shù)r(x,α)；

(4)根據(jù)似然函數(shù)l(x)和懲罰函數(shù)r(x,α)建立pet重建目標函數(shù)ψ(x,α)如下，通過對其進行最優(yōu)化求解重建得到pet濃度分布矩陣x并依此進行pet成像；

ψ(x,α)＝λl(x)+r(x,α)

其中：λ為權重系數(shù)，α為稀疏系數(shù)矩陣。

所述似然函數(shù)l(x)的表達式如下：

其中：yi為符合計數(shù)矩陣y中對應第i組探測器對所采集到的符合計數(shù)向量，為符合計數(shù)向量yi中所有元素的平均值，m為探測器對的數(shù)量。

所述懲罰函數(shù)r(x,α)的表達式如下：

其中：es為分割算子，esx為經(jīng)分割后pet濃度分布矩陣x中的第s個塊矩陣，μ為權重系數(shù)，稀疏系數(shù)矩陣α的大小為m×n，n＝(m-n+1)²，m為pet濃度分布矩陣x的維度，n為預設的塊矩陣維度，αs為稀疏系數(shù)矩陣α中第s列稀疏系數(shù)向量，||αs||0表示稀疏系數(shù)向量αs中非零元素的個數(shù)，||||2表示二范數(shù)。

所述塊字典d由解剖圖像通過塊字典學習算法訓練得到，所述塊字典學習算法包括稀疏凝聚算法求解字典塊結構、塊形式的正交匹配追蹤算法求解稀疏系數(shù)以及塊形式的字典元素更新。

所述塊字典學習算法基于以下函數(shù)求解得到塊字典d：

s.t.||xs||0,d≤k且|dj|≤smax,j∈d

其中：y為用于訓練的解剖圖像且大小為n²×n，塊字典d的大小為n²×l，l為設定值即塊字典d的列數(shù)，d為塊字典d列聚類后得到對應列的類別標簽數(shù)組，x為稀疏系數(shù)矩陣且大小為l×n，xs為稀疏系數(shù)矩陣x中的第s列稀疏系數(shù)向量，||xs||0,d表示稀疏系數(shù)向量xs中連續(xù)非0元素組的個數(shù)，k為設定的稀疏度，|dj|為類別標簽數(shù)組d中類別標簽為j的個數(shù)，j為自然數(shù)且1≤j≤k，k為設定聚類后的類別數(shù)，所述連續(xù)非0元素組中的元素個數(shù)小于等于smax，smax為設定的元素個數(shù)上限值。

優(yōu)選地，所述解剖圖像y選用與pet圖像對象具有相似性的ct圖像。

本發(fā)明通過塊字典學習的方法從解剖圖像中學習先驗知識，利用稀疏表達建立了基于先驗知識的懲罰項，并結合根據(jù)泊松分布模型建立的似然函數(shù)形成后驗函數(shù)形式的目標方程；其中，本發(fā)明利用解剖圖像和塊字典學習方法訓練的塊字典可以代表解剖圖像的局部特征，并且稀疏度相同的情況下，表現(xiàn)出比未采用塊結構字典更好的性能。因此，本發(fā)明pet重建方法有效地克服了ml-em方法中出現(xiàn)的棋盤效應等病態(tài)解問題，同時相比于其他重建方法，也取得比較好的重建效果。

附圖說明

圖1(a)為本發(fā)明pet圖像重建方法的流程示意圖。

圖1(b)為本發(fā)明塊字典學習方法的流程示意圖。

圖2為蒙特卡洛模擬zubal胸腔數(shù)據(jù)的模板示意圖。

圖3(a)為zubal胸腔數(shù)據(jù)的真值圖像。

圖3(b)為采用ml-em方法對zubal胸腔數(shù)據(jù)的重建結果示意圖。

圖3(c)為采用本發(fā)明方法對zubal胸腔數(shù)據(jù)的重建結果示意圖。

圖3(d)為采用sps-os方法對zubal胸腔數(shù)據(jù)的重建結果示意圖。

圖4為蒙特卡洛模擬探測目標數(shù)據(jù)的模板示意圖。

圖5(a)為探測目標數(shù)據(jù)的真值圖像。

圖5(b)為采用ml-em方法對探測目標數(shù)據(jù)的重建結果示意圖。

圖5(c)為采用sps-os方法對探測目標數(shù)據(jù)的重建結果示意圖。

圖5(d)為采用本發(fā)明方法對探測目標數(shù)據(jù)的重建結果示意圖。

圖6(a)為對探測目標數(shù)據(jù)真值圖像的聚類結果示意圖。

圖6(b)為對采用ml-em方法重建圖像的聚類結果示意圖。

圖6(c)為對采用sps-os方法重建圖像的聚類結果示意圖。

圖6(d)為對采用本發(fā)明方法重建圖像的聚類結果示意圖。

具體實施方式

為了更為具體地描述本發(fā)明，下面結合附圖及具體實施方式對本發(fā)明的技術方案進行詳細說明。

如圖1(a)所示，本發(fā)明基于塊字典學習和稀疏表達的pet圖像重建方法，包括如下步驟：

(1)利用pet成像儀器對注入放射性示蹤劑的生物組織進行探測，通過各個位置探測器對收集到符合計數(shù)向量，建立符合計數(shù)矩陣y；

(2)根據(jù)pet成像原理以及pet成像儀器的系統(tǒng)參數(shù)，建立pet測量方程如下：

y＝gx+r+v

其中：g為系統(tǒng)矩陣，由pet測量系統(tǒng)決定，其中的元素gij表示的是x中像素j被探測器i探測到的概率，r為因反射符合事件引起的測量噪聲，v為因散射符合事件引起的測量噪聲。

(3)根據(jù)每個探測對上符合計數(shù)事件符合泊松分布的假設，建立基于泊松分布的似然函數(shù)l(x)如下：

其中：y＝{yi,i＝1,2…m}，是系統(tǒng)每個探測器對探測到的符合事件的個數(shù)和，m是探測器對的個數(shù)。每個探測器探測到的光子個數(shù)都是獨立同分布，且符合伯努利過程，正弦圖中每個元素都可以用一個符合泊松分布的隨機變量來描述如下：

(4)用塊字典學習方法根據(jù)解剖圖像訓練塊字典，并根據(jù)其稀疏表達，形成懲罰函數(shù)r(x，α)如下：

其中：αs為第s塊區(qū)域的稀疏表達系數(shù)，μ為稀疏項的權重系數(shù)。每個字典元素的長度為n²，將圖像分成n×n的小塊，并且用一個向量表示。用一個分割算子e從圖像向量中提取區(qū)域塊。將圖像排列成的向量，并且將其分成n個重疊的區(qū)塊，則代表每一個小塊向量可以表示成ps＝esx，d為由解剖圖像通過塊字典學習方法訓練得到的塊字典。

塊字典d由解剖圖像用塊字典學習算法得到，如圖1(b)所示，塊字典學習算法包括稀疏凝聚算法求解字典塊結構，塊形式的正交匹配追蹤算法求解稀疏系數(shù)，塊形式的字典元素更新方法；塊字典學習算法主要基于以下函數(shù)表達：

s.t.||xi||0,d≤k,1≤i≤n,|dj|≤smax,j∈d

其中：y為輸入的訓練信號，x為稀疏系數(shù)矩陣，xi表示稀疏系數(shù)矩陣x的第i列，k表示塊稀疏度。smax為字典中塊的大小的上限，即塊中最多含有的字典元素，|dj|表示字典中第j塊中包含的字典原子個數(shù)，即塊的大小。

塊結構dj由稀疏凝聚算法(sac)得到，稀疏系數(shù)矩陣x由塊正交匹配追蹤算法(bomp)得到，然后以塊形式由k次奇異值分解方法更新字典元素，其中稀疏凝聚算法將對字典元素的分組轉換為對相應的系數(shù)矩陣的行的分組，將具有相似非零元素分布的行歸為一組等效于將字典中相似的列歸為一塊：

其中，wj(x,d)表示系數(shù)矩陣x中塊dj中的字典元素所對應的行是非零值的列集；設置最大塊尺寸s,求解以下式子：

其中，|wj|表示的wj的大小，b是指字典當前的塊數(shù)量；在滿足|dj|≤s的前提下，合并相似的塊，每一次迭代中，找到一對塊滿足：

合并和即將wj1←{wj1∪wj2}和重復這一步，直到?jīng)]有塊可以合并。

然后，利用塊形式的正交匹配追蹤算法更新稀疏系數(shù)矩陣，具體地：設置塊稀疏度和容忍誤差，輸入字典d和訓練信號矩陣y，按照以下所示步驟更新系數(shù)；塊稀疏度和容忍誤差都需要根據(jù)具體的圖像而定，塊稀疏度一般設置為2或3，而容忍誤差一般為圖像中像素均值的0.1左右。

其中：r為信號y和稀疏表達的殘差；上述算法得到的稀疏系數(shù)矩陣x以及塊結構d，結合原字典d，對于字典中的某塊dj，將誤差函數(shù)寫成然后用svd奇異值分解這里取秩為|dj|的對角矩陣即得到塊dj更新后的字典元素和相應系數(shù)如下；將字典中所有的塊更新之后，即得到新的字典d。

(5)將似然函數(shù)和懲罰函數(shù)相結合，形成新的目標函數(shù)如下：

ψ(x,α)＝λl(x)+r(x,α)

其中：λ為懲罰函數(shù)和似然函數(shù)之間的關系系數(shù)，α為稀疏系數(shù)矩陣。

(6)對目標函數(shù)ψ(x,α)最優(yōu)化求解即得到pet濃度分布矩陣x，目標函數(shù)ψ(x,α)的優(yōu)化求解表達式如下：

針對上式，首先用期望值最大化的思想，引入一個隱含變量cij，表示像素j對于探測器i貢獻的符合事件，gij為系統(tǒng)矩陣的元素，而xj表示像素j發(fā)出的光子對數(shù)，進一步表達如下：

計算出cij的期望值再將其帶回目標方程得到新的目標方程

將新的目標函數(shù)的第二項用可分離的凸替代函數(shù)替換：

進一步將[esx]l分解成含xj的形式：

其中，為第k次迭代中對重建圖像第j個元素的估計值，es,lj為es,l的第j個元素；利用以下不等式，即得到替代的凸函數(shù)。

新的目標方程變成：

在以上的目標方程中，對xj求偏導，即得到的交替計算em中的期望值和最值直至收斂。

以下我們通過對蒙特卡洛模擬數(shù)據(jù)進行實驗來驗證本發(fā)明pet重建的準確性。實驗運行環(huán)境為:8g內(nèi)存，2.70ghz，64位操作系統(tǒng)，cpu為intelcorei5；所模擬的pet掃描儀型號為hamamatsushr-22000，設定的放射性核素及藥物為18f-fdg，設置正弦圖為128個投影角度在每個角度下128條射束采集到的數(shù)據(jù)結果，系統(tǒng)矩陣g∈r^16384×16384；實驗中使1×10⁶計數(shù)率下的投影數(shù)據(jù)作為實驗數(shù)據(jù)。

將本發(fā)明pet重建方法和經(jīng)典的ml-em方法以及sps-os方法進行了比較，為控制結果的可比性，使用了相同的測量數(shù)據(jù)y和系統(tǒng)矩陣g。圖2為實驗用的zubal胸腔數(shù)據(jù)的模板，除背景外將圖像分為三個感興趣的區(qū)域roi；圖3(a)為真值圖，圖3(b)～圖3(d)分別為ml-em，本發(fā)明方法以及sps-os方法的重建效果。計算每個roi區(qū)域相當于真值圖像的偏差和方差，表1是分析結果：

表1

從結果來看，本發(fā)明重建框架不僅比傳統(tǒng)ml-em的方法重建精度更高，更平滑；相比于sps-os方法，本發(fā)明重建框架的重建效果也有較好的表現(xiàn)。

圖4為實驗用的6個尺寸不同的探測目標的模板，圖5(a)為真值圖，圖5(b)～圖5(d)分別為ml-em，sps-os和本發(fā)明方法的重建效果圖；圖6(a)～圖6(d)是相對應的經(jīng)過k均值聚類后的結果，用jaccard指數(shù)來分析每一個探測目標的探測效果，如表2所示：

表2

從結果來看，本發(fā)明重建框架的探測性能與sps-os的探測性能相當，均優(yōu)于傳統(tǒng)的ml-em方法。雖然傳統(tǒng)ml-em方法能探測到最小的目標，但是從重建的圖中可以看到，目標周圍有許多噪聲，不適合實際中使用。

上述對實施例的描述是為便于本技術領域的普通技術人員能理解和應用本發(fā)明。熟悉本領域技術的人員顯然可以容易地對上述實施例做出各種修改，并把在此說明的一般原理應用到其他實施例中而不必經(jīng)過創(chuàng)造性的勞動。因此，本發(fā)明不限于上述實施例，本領域技術人員根據(jù)本發(fā)明的揭示，對于本發(fā)明做出的改進和修改都應該在本發(fā)明的保護范圍之內(nèi)。