本發(fā)明屬于pet成像技術領域，具體涉及一種基于圖像分割的動態(tài)pet圖像重建及示蹤動力學參數(shù)估計方法。

背景技術：

正電子發(fā)射斷層成像(positronemissiontomography，簡稱pet)是核醫(yī)學成像的一種，在生物醫(yī)學研究和臨床診療中正逐漸顯示出其重要的作用。pet成像技術通過對注入生物體內(nèi)示蹤劑的放射性分布進行顯像，能夠從分子水平上發(fā)現(xiàn)細胞的代謝異常，為疾病的早期診斷和預防提供了有效依據(jù)。相對于只能提供一定時間窗內(nèi)生物體組織中放射性強度的平均分布的靜態(tài)pet成像技術而言，動態(tài)pet成像基于動力學模型能夠獲得描述生物體生命活動的定量功能參數(shù)，在科學研究和臨床應用中具有重要應用價值。

在動態(tài)pet的重建問題中，通常的做法是將圖像分割、圖像重建與動力學參數(shù)估計作為三個單獨的問題分開求解。隨著動態(tài)pet成像技術對時間分辨率的要求越來越高，每一幀測量數(shù)據(jù)中有限的光子計數(shù)不斷地挑戰(zhàn)著動態(tài)pet成像質量。此時，由于單幀數(shù)據(jù)的光子計數(shù)不足無法較好地反映出統(tǒng)計特性，利用基于統(tǒng)計特性的重建算法對每一幀數(shù)據(jù)進行單獨重建的方法往往無法得到十分準確的重建結果。此時，通過結合示蹤劑動力學模型引入來自時間縱軸的放射性強度分布的先驗信息，能夠在提高pet重建圖像質量的基礎上同時實現(xiàn)了對動力學參數(shù)的估計。此外，圖像分割也是pet成像技術中的一個傳統(tǒng)難題。除了對pet圖像進行分割以外，還有一類做法是對動態(tài)pet的時間放射性曲線(time-activitycurve，tac)來進行聚類的操作，從而區(qū)分出生物組織里不同的功能區(qū)域。但無論是從pet的圖像平面上還是從動力學參數(shù)的角度來看待這個問題，基于重建結果的分割算法準確度始終依賴與前一步重建結果的精確度，無法解決分割算法對測量噪聲敏感的問題。

實際上，圖像分割、圖像重建與動力學參數(shù)估計這三個問題之間的關系是十分緊密的。通過將三者耦合進入同一個目標方程，通過選擇合適的噪聲模型能夠更好地匹配pet測量數(shù)據(jù)，此時分割結果對噪聲的敏感度將會降低；而分割結果的引入又能夠使已有的圖像估計結果更加準確，從而總體上獲得一個更加符合真實情況的重建結果。

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

鑒于上述，本發(fā)明提供了一種基于圖像分割的動態(tài)pet圖像重建及示蹤動力學參數(shù)估計方法，能夠同時實現(xiàn)功能區(qū)域的分割以及動態(tài)pet聯(lián)合重建。

一種基于圖像分割的動態(tài)pet圖像重建及示蹤動力學參數(shù)估計方法，包括如下步驟：

(1)利用探測器對注入有放射性藥劑的生物組織進行探測，動態(tài)采集得到對應各個時刻的符合計數(shù)向量，并建立符合計數(shù)矩陣y；

(2)使動態(tài)的pet圖像序列組合成pet濃度分布矩陣x，根據(jù)pet成像原理，建立pet測量方程；

(3)通過對pet測量方程引入全變分(totalvariation，tv)約束，得到基于tv的pet圖像重建模型l(x)；

(4)利用房室模型匹配估計示蹤動力學參數(shù)，建立關于x和ф的同步重建模型s(x，ф)；

(5)基于預處理得到的區(qū)域分割結果，對示蹤動力學參數(shù)矩陣ф進行聚類，得到聚類分割模型c(ф)；

(6)將上述三個模型l(x)、s(x，ф)和c(ф)相結合得到同步重建的目標函數(shù)lsc(x，ф)如下：

其中：γ和∈均為權重系數(shù)；

(7)對目標函數(shù)lsc(x，ф)進行最優(yōu)化求解后即得到pet濃度分布矩陣x和示蹤動力學參數(shù)矩陣ф。

所述符合計數(shù)矩陣y由各符合計數(shù)向量按時序排列組成，所述pet濃度分布矩陣x由各時刻對應的pet濃度分布向量(即一幀pet圖像)按時序排列組成。

所述pet測量方程的表達式如下：

y＝gx+r+s

其中：g為系統(tǒng)矩陣，r和s分別為反映隨機事件和散射事件的測量噪聲矩陣。

所述pet圖像重建模型l(x)的表達式如下：

其中：α為權重系數(shù)，tv(x)為關于x的全變分正則項，gij為系統(tǒng)矩陣g中第i行第j列元素值(其表示的是第j個像素點處出射的光子被第i個探測器接收到的概率)，yim為符合計數(shù)矩陣y中第i行第m列元素值，xjm為pet濃度分布矩陣x中第j行第m列元素值，rim為反映隨機事件的測量噪聲矩陣r中第i行第m列元素值，sim為反映散射事件的測量噪聲矩陣s中第i行第m列元素值，i、j和m均為自然數(shù)且1≤i≤n，1≤j≤k，1≤m≤m，n為符合計數(shù)向量的維度，k為pet濃度分布矩陣x的行數(shù)即pet圖像的像素點個數(shù)，m為pet濃度分布矩陣x的列數(shù)即采樣時間長度。

所述全變分正則項tv(x)的表達式如下：

其中：d(xjm)為關于xjm的二維差分向量，該向量第一行元素值為xjm-xj，m+1，第二行元素值為xjm-xj+1，m，xj，m+1為pet濃度分布矩陣x中第j行第m+1列元素值，xj+1，m為pet濃度分布矩陣x中第j+1行第m列元素值，||||2表示2范數(shù)。

所述同步重建模型s(x，ф)的表達式如下：

其中：μ為權重系數(shù)，ψ為字典矩陣，^t表示轉置，||||2表示2范數(shù)，||||1表示1范數(shù)。

所述聚類分割模型c(ф)的表達式如下：

其中：ch為示蹤動力學參數(shù)矩陣ф中屬于第h類的參數(shù)向量集合，φh為參數(shù)向量集合ch中的任一參數(shù)向量，nh為參數(shù)向量集合ch中的參數(shù)向量個數(shù)，h為自然數(shù)且1≤h≤h，h為聚類的類數(shù)，^t表示轉置。

所述步驟(7)中采用admm(alternatingdirectionmethodofmultipliers，交替方向乘子算法)結合軟閾值(soft-thresholding)迭代優(yōu)化算法對目標函數(shù)lsc(x，ф)進行最優(yōu)化求解；其中，admm針對pet濃度分布矩陣x進行迭代優(yōu)化求解，軟閾值針對示蹤動力學參數(shù)矩陣ф迭代優(yōu)化求解。

所述軟閾值迭代優(yōu)化算法基于以下方程進行線性處理以對示蹤動力學參數(shù)矩陣ф進行迭代更新：

nh＝ih-eh/nh

其中：eh為一個全部由1組成的nh×nh維矩陣，ih為一個與矩陣eh維度相同的單位矩陣，μ為權重系數(shù)，^t表示轉置，ψ為字典矩陣，ch為示蹤動力學參數(shù)矩陣ф中屬于第h類的參數(shù)向量集合，φh為參數(shù)向量集合ch中的任一參數(shù)向量，nh為參數(shù)向量集合ch中的參數(shù)向量個數(shù)，xh為pet濃度分布矩陣x中與參數(shù)向量φh對應的一條tac，||||2表示2范數(shù)，||||1表示1范數(shù)，h為自然數(shù)且1≤h≤h，h為聚類的類數(shù)。

所述字典矩陣ψ的表達式如下：

其中：ci(t)和ci(τ)分別為t時刻和τ時刻放射性藥劑在血漿中的濃度值，和分別為關于第m組符合計數(shù)向量采集的開始時間和結束時間，θc對應為第c個房室組織指數(shù)函數(shù)的系數(shù)，m和c均為自然數(shù)且1≤m≤m，1≤c≤z，m為pet濃度分布矩陣x的列數(shù)即采樣時間長度，z為大于1的自然數(shù)；θ1～θn的取值為在區(qū)間[θmin，θmax]內(nèi)按指數(shù)級間隔進行選取，θmin和θmax分別為系數(shù)的上下限閾值，t和τ均表示時間。

本發(fā)明通過將統(tǒng)計學重建模型與有生理意義的模型相結合能夠得到更加準確的聯(lián)合重建結果。當重建和分割耦合到一個聯(lián)合或同時的求解框架中時，對分割任務而言，可以獲得噪聲模型對原始投影數(shù)據(jù)進行建模的信息，而基于分割的結果也能夠增強每個區(qū)域內(nèi)的均勻性，從而實現(xiàn)一個更符合真實情況的重建結果。與其他單獨重建動態(tài)pet圖像或估計動力學參數(shù)的算法相比，本發(fā)明也都能獲得較好的重建結果。結合本發(fā)明在模擬數(shù)據(jù)和真實數(shù)據(jù)實驗中的表現(xiàn)，與其他單獨重建動態(tài)pet圖像、估計動力學參數(shù)以及圖像分割的算法相比，本發(fā)明也都能獲得較好的重建結果。

附圖說明

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

圖2(a)為zubal胸腔數(shù)據(jù)第4幀的真實圖像。

圖2(b)為數(shù)據(jù)計數(shù)率為5*10^4下采用ml-em方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第4幀圖像結果。

圖2(c)為數(shù)據(jù)計數(shù)率為1*10^5下采用ml-em方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第4幀圖像結果。

圖2(d)為數(shù)據(jù)計數(shù)率為5*10^4下采用本發(fā)明方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第4幀圖像結果。

圖2(e)為數(shù)據(jù)計數(shù)率為1*10^5下采用本發(fā)明方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第4幀圖像結果。

圖3(a)為zubal胸腔數(shù)據(jù)第6幀的真實圖像。

圖3(b)為數(shù)據(jù)計數(shù)率為5*10^4下采用ml-em方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第6幀圖像結果。

圖3(c)為數(shù)據(jù)計數(shù)率為1*10^5下采用ml-em方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第6幀圖像結果。

圖3(d)為數(shù)據(jù)計數(shù)率為5*10^4下采用本發(fā)明方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第6幀圖像結果。

圖3(e)為數(shù)據(jù)計數(shù)率為1*10^5下采用本發(fā)明方法對蒙特卡洛模擬zubal胸腔數(shù)據(jù)重建的第6幀圖像結果。

圖4(a)為蒙特卡洛模擬zubal胸腔數(shù)據(jù)在采用本發(fā)明方法下的聚類結果示意圖。

圖4(b)為蒙特卡洛模擬zubal胸腔數(shù)據(jù)在采用k均值分類方法下的聚類結果示意圖。

圖4(c)為蒙特卡洛模擬zubal胸腔數(shù)據(jù)在采用基于動力學模型譜分類方法下的聚類結果示意圖。

具體實施方式

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

本發(fā)明基于區(qū)域分割的動態(tài)pet及動力學參數(shù)聯(lián)合重建方法，包括如下步驟：

(1)根據(jù)動態(tài)pet掃描的模式建立測量數(shù)據(jù)矩陣y和系統(tǒng)矩陣g。

將動態(tài)pet的掃描過程按照需要劃分出一定數(shù)量的時間幀，每個時間幀內(nèi)探測器采集得到的符合計數(shù)向量的可以按照時間的順序構建出一個動態(tài)pet的測量數(shù)據(jù)矩陣y；并統(tǒng)計每一像素點處出射的光子被各探測器接收到的概率，從而得到系統(tǒng)矩陣g。

(2)根據(jù)測量數(shù)據(jù)的標記示蹤化合物的正電子核素種類以及時間幀間隔的分布，基于雙房室模型理論建立時間基函數(shù)的字典矩陣ψ。

若是模擬數(shù)據(jù)，則根據(jù)模擬的示蹤同位核素的半衰期設置指數(shù)參數(shù)θ和探測器掃描時間間隔建立相應的字典矩陣；若是真實數(shù)據(jù)，則按照實際操作中使用的示蹤同位核素和實驗設置建立相應的字典矩陣。

根據(jù)房室模型可以解出組織中的放射性強度隨時間的變化，即時間-放射性函數(shù)是輸入組織的tac和一個δ(t)函數(shù)以及一組指數(shù)函數(shù)的卷積，即是所有組織房室的tac之和。因此，一個pet示蹤房室模型可以用一組一階線性方程來表示，每個方程表示一個不同的房室中的示蹤劑放射性強度，即：

s.t.ψ0(t)＝ci(t)

其中：ct(t)表示t時刻組織中的放射性強度。n表示不同房室的總數(shù)，ψc(t)表示第c種房室對應的放射性強度隨時間變化的函數(shù)，φc是ψc(t)相對應的系數(shù)。θc是對應房室的指數(shù)函數(shù)中的參數(shù)。ci(t)表示t時刻輸入的放射性強度。

重建的每一幀的動態(tài)pet圖像可看作在這一幀與上一幀的時間間隔內(nèi)對房室組織時間-放射性曲線的積分，即：

其中，xjm是pet圖像序列中第m幀的像素j所對應的放射性強度，和分別是第m幀的開始時間和結束時間，是方式組織的像素j處對應的tac。基于此關系，我們將一組tac對應的方程組擴展為一個字典矩陣如下：

其中：

(3)初始化，設置權重系數(shù)α、γ和μ，軟閾值算法中的步長σ，迭代次數(shù)k置0，設置最大迭代次數(shù)kmax。

權重系數(shù)α在[2^7.5,2^8.5]左右，這個是根據(jù)原來已有的算法給出的參考范圍進行調(diào)整得到的。權重系數(shù)γ一般在[1,5]左右的范圍內(nèi)，這個是基于實驗的結果獲得的。權重系數(shù)μ一般可在[0.01,0.15]的范圍內(nèi)進行選擇得到合適的值(根據(jù)參考文獻positronemissiontomographycompartmentalmodels&colon；abasispursuitstrategyforkineticmodeling中的實驗結果)，當噪聲增加的時候，算法不能準確地從字典中確定出正確的時間放射性曲線，此時需要適當增大μ的取值，增強稀疏約束。

(4)進入算法迭代更新過程，對于第k次迭代,固定動力學參數(shù)φ^k得到pet圖像的更新結果x^k+1。

4.1提取出目標參數(shù)中僅包括pet圖像x的部分獲得以下子優(yōu)化問題：

其中：tv(x)為全變分正則項，α為權重系數(shù)，i＝1，…，n為探測器的編號，m＝1，…，m代表時間幀，j＝1，…，k代表成像平面上像素點的編號；因此yim是第m幀第i個探測器的測量數(shù)據(jù)，xjm是第m幀第j個像素點的濃度值，而gij則是系統(tǒng)矩陣的一個元素，代表第j個像素點處出射的光子被第i個探測器接收到的概率，rim和sim分別表示相對應的時間幀和探測器測量數(shù)據(jù)中的隨機事件和散射事件。

4.2基于全變分正則項的計算公式定義新變量ωjm并將4.1中的目標方程轉化為增廣拉格朗日方程如下：

其中：其中向量υjm是拉格朗日乘子，βjm則是懲罰系數(shù)，xm是第m幀pet圖像像素點組成的向量，則是xm在第j個像素點處的離散偏微算子。

4.3固定x，利用二維收縮公式更新變量ωjm；

4.4固定ωjm，將la(ωjm，x)對x求導后按照如下公式更新pet圖像x：

x^k+1＝x^k-ρ^kd^k(x)x^k

其中：ρ^k是通過逆向非單調(diào)線搜索確定的最快下降步長，d^k(x)是增廣拉格朗日方程對x的導數(shù)；

4.5固定x和ωjm，更新拉格朗日乘子；

4.6判斷是否滿足迭代停止條件：x的變化＜10^-3或達到最大子迭代次數(shù)，不滿足則返回4.4；若滿足則迭代停止，得到pet圖像x的更新結果x^k+1。

(5)對于第k次迭代,固定pet圖像x^k+1和動力學參數(shù)φ^k，使用狄利克雷聚類過程更新聚類結果

(6)對于第k次迭代,固定pet圖像x^k+1得到動力學參數(shù)的估計結果φ^k+1：

6.1提取出目標參數(shù)中僅包括動力學參數(shù)φ的部分獲得以下子優(yōu)化問題：

其中：代表著系數(shù)矩陣φ中屬于h類的像素點的子集，nh是該集合中像素點的個數(shù)，ρh由第h類對應的像素點的平均值組成的矩陣，∈為控制子模型的權重系數(shù)。

6.2通過定義一個新的矩陣nh＝ih-eh/nh，假設eh是一個全部由1組成的nh×nh維的矩陣，ih是一個與eh維度相同的單位矩陣并定義新矩陣nh＝ih-eh/nh。引入化簡上述方程：

6.3對上式進行線性化處理后利用軟閾值迭代算子可以得到更新的動力學參數(shù)φ^k+1。

(7)判斷是否滿足迭代停止條件：x和φ的變化＜10^-3或達到最大迭代次數(shù)kmax，不滿足則返回步驟(4)，若滿足則迭代停止，得到pet圖像x、示蹤動力學參數(shù)φ以及聚類結果。

以下我們通過對蒙特卡洛模擬的zubal胸腔模板數(shù)據(jù)進行實驗從而驗證本發(fā)明系統(tǒng)重建和分割結果的準確性，圖1為實驗所用的zubal胸腔數(shù)據(jù)的模板示意圖，將不同的區(qū)域分為三個感興趣的區(qū)域(regionofinterest，roi)。實驗運行環(huán)境為：8g內(nèi)存，3.40ghz，64位操作系統(tǒng)，cpu為inteli7-3770；所模擬的pet掃描儀型號為hamamatsushr-22000，設定的放射性核素及藥物為¹⁸f-fdg，設置sinogram為128個投影角度在每個角度下128條射束采集到的數(shù)據(jù)結果，系統(tǒng)矩陣g的大小為16384×16384。在本次試驗中，對5×10⁴、1×10⁵兩種不同的計數(shù)率下的投影數(shù)據(jù)進行實驗。

本發(fā)明重建框架中對pet圖像的結果和傳統(tǒng)的ml-em(最大似然-期望最大法)的重建結果做比較，二者使用相同的測量數(shù)據(jù)矩陣y和系統(tǒng)矩陣g以控制結果的可比性。圖2(a)和圖3(a)均為分布真值，圖2(b)～圖2(e)分別是在ml-em在計數(shù)率為5×10⁴和1×10⁵的情況下以及本發(fā)明重建框架在計數(shù)率為5×10⁴和1×10⁵的情況下的動態(tài)pet圖像序列中第4幀的重建結果；圖3(b)～圖3(e)分別是在ml-em在計數(shù)率為5×10⁴和1×10⁵的情況下以及本發(fā)明重建框架在計數(shù)率為5×10⁴和1×10⁵的情況下的動態(tài)pet圖像序列中第6幀的重建結果。以明顯看出聯(lián)合重接獲得的結果在區(qū)域內(nèi)噪聲明顯小于ml-em獲得的圖像，在保證邊緣對比的情況下功能區(qū)域內(nèi)更加平滑，表1是其進一步的量化分析結果。

表1

圖4(a)～圖4(c)是聚類的比較結果，本發(fā)明重建框架得到的功能區(qū)域分布和已有的k均值分類以及基于動力學模型的譜分類方法得到的結果進行比較?？梢钥闯雎?lián)合重建的聚類結果中的噪聲明顯較大，這主要是由于狄利克雷過程在進行分類預處理時是以像素點為單位進行處理的，因此對噪聲十分敏感。

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