一種基於圖像分割的動態PET圖像重建及示蹤動力學參數估計方法與流程

2023-10-11 15:08:44 2

本發明屬於pet成像技術領域，具體涉及一種基於圖像分割的動態pet圖像重建及示蹤動力學參數估計方法。

背景技術：

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

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

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

技術實現要素：

鑑於上述，本發明提供了一種基於圖像分割的動態pet圖像重建及示蹤動力學參數估計方法，能夠同時實現功能區域的分割以及動態pet聯合重建。

一種基於圖像分割的動態pet圖像重建及示蹤動力學參數估計方法，包括如下步驟：

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

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

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

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

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

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

其中：γ和∈均為權重係數；

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

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

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

y＝gx+r+s

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

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

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

所述全變分正則項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範數。

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

其中：μ為權重係數，ψ為字典矩陣，t表示轉置，||||2表示2範數，||||1表示1範數。

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

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

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

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

nh＝ih-eh/nh

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

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

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

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

附圖說明

圖1為蒙特卡洛模擬zubal胸腔數據的模板圖像。

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

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

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

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

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

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

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

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

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

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

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

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

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

具體實施方式

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

本發明基於區域分割的動態pet及動力學參數聯合重建方法，包括如下步驟：

(1)根據動態pet掃描的模式建立測量數據矩陣y和系統矩陣g。

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

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

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

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

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

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

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

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

其中：

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

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

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

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

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

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

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

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

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

xk+1＝xk-ρkdk(x)xk

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

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

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

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

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

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

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

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

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

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

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

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

表1

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

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