本發(fā)明屬于醫(yī)學(xué)影像領(lǐng)域，涉及統(tǒng)計輪廓模型和圖割分割。

背景技術(shù)：

醫(yī)學(xué)圖像分割作為一種可以在非侵入的情況下提供人體結(jié)構(gòu)信息的技術(shù)，是醫(yī)學(xué)影像數(shù)據(jù)分析和可視化的第一步，被廣泛用于可視化學(xué)習(xí)人體結(jié)構(gòu)的解剖學(xué)特性、模擬生理過程、確定病理位置、跟蹤疾病發(fā)展過程、以及確定放射治療劑量或手術(shù)規(guī)模。因此，醫(yī)學(xué)圖像分割是計算機輔助治療系統(tǒng)中關(guān)鍵的一步，分割的準確性直接影響輔助治療的效果。

但是，醫(yī)學(xué)圖像的準確分割面臨許多挑戰(zhàn)。首先，在像素或體素灰度的空間重疊下，許多解剖學(xué)結(jié)構(gòu)是非均一的。其次，醫(yī)學(xué)圖像中存在低對比度部分，比如腎臟和心臟的邊緣在CT圖像中都較難辨別，圖像采集過程中的噪音會進一步降低低對比度部分的辨識度。最后，當(dāng)采用三維圖像時，圖像采集過程中目標器官會產(chǎn)生形狀的內(nèi)部和主體變異，這使得同一器官在圖像中展示出的形態(tài)千差萬別。

為了克服這些困難，醫(yī)學(xué)圖像分割采用了多種分割方法，但總體上，沒有一個分割方法可以對任何器官都獲得較好的分割結(jié)果。最基本的閾值法簡單快速，但對與背景灰度差異不明顯的器官無法較好的分割，且對噪聲十分敏感。區(qū)域增長法和水閥法和基于閾值的方法有相同的缺陷，更容易出現(xiàn)過分割的情況。近年來，基于圖割(Graph Cut，簡稱GC)的方法在醫(yī)學(xué)圖像分割中展示出較好的分割結(jié)果，和主動輪廓、水平集方法類似，GC是基于能量的方法：將一副圖像轉(zhuǎn)化成一個圖(Graph)，進而把圖像分割問題轉(zhuǎn)化為能量方程的最小化問題。GC方法的優(yōu)勢在于把圖像繪圖和邊界信息融合到能量方程中，認為當(dāng)一個圖像得到恰當(dāng)分割時，能量方程的值最小。

上述方法都常被直接用于圖像數(shù)據(jù)，沒有考慮待分割目標的形狀先驗知識，而這些形狀先驗知識在醫(yī)學(xué)圖像的器官分割中是普遍存在的。分割問題中，可以通過兩種方法表示形狀先驗信息：主動形狀模型(active shape model,ASM)和主動外觀模型(active appearance model，AAM)。把手動分割的圖像作為訓(xùn)練集，對其進行主成分分析來獲取平均形狀和形狀變化，這樣訓(xùn)練集中的所有形狀都可以用平均形狀和形狀變化的線性組合來表示。雖然主動形狀模型或主動外觀模型可以得到可接受的分割結(jié)果，但這個分割結(jié)果往往不是全局最優(yōu)的。此外，單純使用形狀信息進行分割，需要足以涵蓋各種形狀變化的訓(xùn)練集，對訓(xùn)練圖像數(shù)目和質(zhì)量有很高要求。

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

本發(fā)明的目的是針對上述醫(yī)學(xué)圖像器官分割存在的分割不準確、效率不高的問題，提出了一種基于統(tǒng)計形狀模型和Graph Cut方法的分割方法。針對Graph Cut分割方法在目標器官與背景差異較小的情況下無法準確分割的問題，本發(fā)明通過引入統(tǒng)計形狀模型加入形狀先驗，結(jié)合待分割器官的統(tǒng)計形狀信息輔助分割。

為了實現(xiàn)上述目的，本發(fā)明包括如下步驟：

(1)建立低對比度器官的統(tǒng)計形狀模型并采集灰度信息，其包括：將包含目標低對比度器官的三維CT圖像作為訓(xùn)練集，手動分割動物體外輪廓、低對比度器官，對其進行有限元離散剖分得到點云數(shù)據(jù)，使用相似性變換將訓(xùn)練集中的樣本對齊配準到一個坐標系上，對配準后的形狀分別使用PCA計算低對比度器官和動物體外輪廓的平均形狀和形狀變化，采集每個特征點的局部灰度信息；并使用機器學(xué)習(xí)方法建立分類函數(shù)，學(xué)習(xí)手動分割的低對比度器官和背景的灰度信息和位置信息；

(2)預(yù)分割低對比度器官，其包括：將步驟(1)中獲得的動物體外輪廓平均形狀的中心與待分割圖像的中心重合，沿特征點法線方向搜索分割得到動物體外輪廓，計算平均形狀到外輪廓的變換T；T可近似為低對比度器官到測試圖像器官的變換，對低對比度器官平均形狀做變換T，得到低對比度器官初始位置T^O；在初始位置T^O基礎(chǔ)上，沿特征點法線方向進行器官的搜索分割，得到低對比度器官預(yù)分割結(jié)果；

(3)初始化Graph，其包括：根據(jù)待分割圖像中像素點的灰度信息和距離信息，以及步驟(1)中獲得的低對比度器官的形狀信息(平均形狀)和分類函數(shù)，構(gòu)造待分割圖像對應(yīng)的Graph；然后使用步驟(2)得到的預(yù)分割結(jié)果初始化Graph中邊的權(quán)值，得到初始化后的Graph，即G_l；

(4)分割低對比度器官，其包括：使用Graph Cut算法對計算步驟(3)中得到的G_l的最大流，即最小化待分割圖像對應(yīng)的能量方程E(L)＝λ(R(L)k₁S_p(L))+(B(L)+k₂S_p1,p2)，得到分割G_l并使得能量方程取值最小的一個割(Cut),即得到最終的低對比度器官邊界。

在上述技術(shù)方案的基礎(chǔ)上，所述步驟(1)中建立動物體外輪廓同低對比度器官統(tǒng)計形狀模型并采集灰度信息的具體步驟為：

(1a)對訓(xùn)練樣本中的低對比度器官和動物體外輪廓進行手動分割，分別對分割結(jié)果進行有限元離散剖分，得到相應(yīng)的三維點云數(shù)據(jù)；對各個樣本的低對比度器官和動物體外輪廓的點云數(shù)據(jù)通過仿射變換進行配準，隨機選取一個樣本的點云數(shù)據(jù)作為基準模板將所有樣本進行旋轉(zhuǎn)、平移等相似性變換，得到兩者對齊配準結(jié)果；

(1b)對步驟(1a)中配準后的的結(jié)果求和平均，分別計算低對比度器官均值模型和動物體外輪廓均值模型

其中，為第i個低對比度器官對齊配準后的樣本，為第i個動物體外輪廓對齊配準后的樣本，L為訓(xùn)練樣本數(shù),均為K行乘3列的矩陣，K為配準后單個樣本三維點云中點的個數(shù)；

(1c)計算低對比度器官配準結(jié)果的協(xié)方差矩陣S^O：

S^O為3K行乘3K列的矩陣；

(1d)對協(xié)方差矩陣S^O進行特征分解，求解對應(yīng)的特征值和特征向量，將特征值按絕對值降序排列，選取前t(1≤t≤3K)個特征值，以確保能夠反應(yīng)樣本變化的主要模式，將主成分對應(yīng)的特征向量單位化處理；

選擇t個特征值使得

其中s_i為第i個奇異值，即選取可以保留原圖像95％的特征；

(1e)根據(jù)步驟(1d)所得，低對比度器官統(tǒng)計形狀模型中的任意形狀X^O均可描述為其均值模型與特征向量P^O和形變參數(shù)b^O的線性組合，即為

為了使形狀變化限制在合理的范圍內(nèi)，形變參數(shù)b^O需滿足：

其中，λ為第一主成分特征值。

(1f)同步驟(1c)-(1e),計算得到動物體外輪廓統(tǒng)計形狀模型中的任意形狀X^C可描述為其均值模型與特征向量P^C和形變參數(shù)b^C的線性組合:

形變參數(shù)b^C需滿足：

其中，λ^C為第一主成分特征值。

(1g)分別提取訓(xùn)練集中動物體外輪廓和低對比度器官所有標記點的局部灰度信息，對每一個標記點沿其輪廓的法線方向在改點兩側(cè)各取n個點，則一個樣本中第i張切片圖像上第j個標記點處的灰度向量為

g_ij＝[g_ij0,g_ij1,...,g_ij2n]^T

再對g_ij進行差分處理以保證偏移的相對不變性及灰度尺度的一致性，差分向量長度為2n，則處理后的灰度向量記為

用上述方法對所有訓(xùn)練樣本中相對應(yīng)的標記點進行處理，得到L個樣本中相對應(yīng)點的平均標準化差分灰度向量，即

計算其協(xié)方差，記為

對所有標記點進行上述操作，得到動物體外輪廓和低對比度器官邊界點的局部灰度信息。

(1h)使用常規(guī)的機器學(xué)習(xí)方法(如Logistic分類)，分別以像素點像素值和坐標為特征，建立兩個分類函數(shù)，學(xué)習(xí)手動分割的低對比度器官和背景的灰度信息和位置信息。

在上述技術(shù)方案的基礎(chǔ)上，所述步驟(2)中通過分割待分割圖像中動物體外輪廓來確定低對比度器官位置，進而使用統(tǒng)計形狀模型預(yù)分割低對比度器官的具體步驟為：

(2a)將步驟(1)中獲得的動物體外輪廓平均形狀的中心與待分割圖像的中心重合，沿特征點法線方向兩側(cè)各k個點共2k+1個點搜索，找到馬氏距離d最小的點即為最佳匹配點，其中

其中q_j為待測邊界點對應(yīng)標記點j的某個特征沿其法線方向的歸一化灰度向量， 和S_yj是標記點j的局部灰度模型。遍歷所有點后產(chǎn)生新的形狀模型，迭代至兩個相鄰的形狀模型差值小于一定閾值時停止搜索。得到動物體外輪廓T^C，計算平均形狀到外輪廓的變換T：

(2b)基于動物體外輪廓與低對比度器官之間的位置相關(guān)性，將配準變換T近似為低對比度器官平均形狀到測試圖像器官的變換，對低對比度器官平均形狀做變換T，得到低對比度器官初始位置P^O，即

(2c)根據(jù)步驟(2b)中確定的低對比度器官初始位置P^O，使用同步驟(2a)所述方法分割得到低對比度器官初始分割結(jié)果T^O；

(2d)對步驟(2c)得到的低對比度器官輪廓T^O的點云數(shù)據(jù)使用三角剖分算法進行面片化，再對得到的三角網(wǎng)格數(shù)據(jù)通過帶標的歐氏距離場進行體素化。

在上述技術(shù)方案的基礎(chǔ)上，所述步驟(3)中根據(jù)待分割圖像構(gòu)造對應(yīng)的Graph,再使用步驟(2)的低對比度器官預(yù)分割結(jié)果初始化Graph的步驟為：

(3a)構(gòu)造Graph：

待分割圖像中的每一個像素點抽象為Graph中的頂點，每個頂點p有兩種邊:

①t-link,邊e_{p,S}，e_{p,T}分別表示像素點和前景、背景的關(guān)系，權(quán)值為

R(L)+k₁S_p(L)

其中k₁為調(diào)節(jié)區(qū)域項R(L)和形狀先驗S_p(L)的正參數(shù)，

R(L)＝-ln Pr(I_r|L)

其中為P_r(I_r|L)為后驗概率，使用步驟(1)中由灰度值得到的分類函數(shù)計算；

Pr_p(L)為點p是被標記為0(背景)或1(低對比度器官)的概率；

當(dāng)時，若點p在內(nèi),否則

當(dāng)時，Pr_p(L)使用步驟(1)由坐標得到的分類函數(shù)計算，T_n為分類函數(shù)決策邊界處像素點坐標與訓(xùn)練樣本分割邊界的平均距離；

②n-link，表示像素點和8鄰居系統(tǒng)中其他八個像素點之間的關(guān)系，權(quán)值為

B(L)+k₂S_p1,p2,

其中k₂為調(diào)節(jié)邊界項B(L)和形狀先驗S_p1,p2的正參數(shù)，

σ為經(jīng)驗值；

其中

L_p1，L_p2為點p₁，p₂的標簽(前景或背景)；為點p_i到平均形狀的歐氏距離；L為訓(xùn)練圖像的數(shù)目，σ_i²由第i個訓(xùn)練圖像計算得到， N_i為第i個訓(xùn)練圖像中像素點的個數(shù)，為點p_j到平均形狀的歐氏距離；

(3b)使用步驟(2)得到的預(yù)分割結(jié)果對步驟(3a)構(gòu)造的Graph進行初始化。根據(jù)步驟(2)預(yù)分割結(jié)果，在低對比度器官內(nèi)部的像素點所對應(yīng)的Graph中頂點的n-link置為e_{p,S}＝Inf,e_{p,T}＝0；在低對比度器官外部的像素點對應(yīng)的n-link置為e_{p,S}＝0,e_{p,T}＝Inf，得到初始化后的圖G_l。

本發(fā)明與現(xiàn)有技術(shù)相比具有如下優(yōu)點：

第一，本發(fā)明使用動物體外輪廓與動物器官的相對位置關(guān)系，利用容易準確分割的動物體外輪廓定位低對比度器官的位置，有效解決統(tǒng)計形狀模型初始位置難以定位的問題。

第二，本發(fā)明結(jié)合目標器官的形狀先驗知識，解決了基于灰度的分割策略難以分割低對比度圖像的問題。

第三，本發(fā)明使用Graph Cut方法，在整幅圖像上計算最小能量值，從而得到全局最優(yōu)解，同時在計算速度上較快。

附圖說明

圖1為本發(fā)明基于統(tǒng)計形狀模型的醫(yī)學(xué)圖像Graph Cut分割方法的基本流程圖；

圖2為本發(fā)明基于統(tǒng)計形狀模型的醫(yī)學(xué)圖像Graph Cut分割方法實施例的流程圖。

具體實施方式

下面結(jié)合附圖詳細說明本發(fā)明技術(shù)方案中所涉及的各個細節(jié)問題。應(yīng)指出的是，所描述的實施例僅旨在便于對本發(fā)明的理解，而對其不起任何限定作用。

在此實施例中，以小鼠腎臟作為分割目標，但不局限于此。實施例的框架如附圖1所示，詳細流程如附圖2所示。

步驟1：獲取小鼠CT斷層數(shù)據(jù)

將已注射造影劑的實驗小鼠固定在Micro-CT成像系統(tǒng)的成像臺上，調(diào)整X射線管、旋轉(zhuǎn)臺以及X射線平板探測器位置，使得三者中心在一條直線上，對小鼠進行360度照射掃描，采集投影數(shù)據(jù)，利用濾波反投影方法對投影數(shù)據(jù)三維重建，得到小鼠CT斷層數(shù)據(jù)。使用3DMed把掃描的CT數(shù)據(jù)進行三維體重建，得到raw格式數(shù)據(jù)。

步驟2：建立腎臟的統(tǒng)計形狀模型并采集灰度信息

(1a)使用MITK Workbench對訓(xùn)練樣本中的腎臟和小鼠體外輪廓進行手動分割，分別對分割結(jié)果進行有限元離散剖分，得到相應(yīng)的三維點云數(shù)據(jù)；對各個樣本的腎臟和小鼠體外輪廓的點云數(shù)據(jù)通過仿射變換進行配準，隨機選取一個樣本的點云數(shù)據(jù)作為基準模板將所有樣本進行旋轉(zhuǎn)、平移等相似性變換，得到兩者對齊配準結(jié)果；

(1b)對步驟(1a)中配準后的的結(jié)果求和平均，分別計算腎臟均值模型和小鼠體外輪廓均值模型

其中，為第i個腎臟對齊配準后的樣本，為第i個小鼠體外輪廓對齊配準后的樣本，L為訓(xùn)練樣本數(shù)，均為K行乘3列的矩陣，K為配準后單個樣本三維點云中點的個數(shù)；

(1c)計算腎臟配準結(jié)果的協(xié)方差矩陣S^O：

S^O為3K行乘3K列的矩陣；

(1d)對協(xié)方差矩陣S^O進行特征分解，求解對應(yīng)的特征值和特征向量，將特征值按絕對值降序排列，選取前t(1≤t≤3K)個特征值，以確保能夠反應(yīng)樣本變化的主要模式，將主成分對應(yīng)的特征向量單位化處理；

選擇t個特征值使得

即選取可以保留原圖像95％的特征；

(1e)根據(jù)步驟(1d)所得，腎臟統(tǒng)計形變模型中的任意形狀X^O均可描述為其均值模型與特征向量P^O和形變參數(shù)b^O的線性組合，即為

為了使形狀變化限制在合理的范圍內(nèi)，形變參數(shù)b^O需滿足：

其中，λ為第一主成分特征值。

(1f)同步驟(1c)-(1e),計算得到小鼠體外輪廓統(tǒng)計形狀模型中的任意形狀X^C可描述為其均值模型與特征向量P^C和形變參數(shù)b^C的線性組合:

形變參數(shù)b^C需滿足：

其中，λ^C為第一主成分特征值。

(1g)分別提取訓(xùn)練集中小鼠體外輪廓和腎臟所有標記點的局部灰度信息，對每一個標記點沿其輪廓的法線方向在改點兩側(cè)各取n個點，則一個樣本中第i張切片圖像上第j個標記點處的灰度向量為

g_ij＝[g_ij0,g_ij1,...,g_ij2n]^T

再對g_ij進行差分處理以保證偏移的相對不變性及灰度尺度的一致性，差分向量長度為2n，則處理后的灰度向量記為

用上述方法對所有訓(xùn)練樣本中相對應(yīng)的標記點進行處理，得到L個樣本中相對應(yīng)點的平均標準化差分灰度向量，即

計算其協(xié)方差，記為

對所有標記點進行上述操作，得到小鼠體外輪廓和腎臟邊界點的局部灰度信息。

(1h)使用機器學(xué)習(xí)方法學(xué)習(xí)手動分割的低對比度器官和背景的灰度信息。本實施例采用Logistic分類，以像素點灰度值和像素點坐標為特征，進行標簽(背景、前景)分類。

步驟3：預(yù)分割腎臟

(2a)將步驟(1)中獲得的小鼠體外輪廓平均形狀的中心與d待分割圖像的中心重合，沿特征點法線方向兩側(cè)各k個點共2k+1個點搜索，找到馬氏距離d最小的點即為最佳匹配點，其中

其中q_j為待測邊界點對應(yīng)標記點j的某個特征沿其法線方向的歸一化灰度向量， 和S_yj是標記點j的局部灰度模型。遍歷所有點后產(chǎn)生新的形狀模型，迭代至兩個相鄰的形狀模型差值小于一定閾值，時停止搜索。得到；小鼠體外輪廓T^C，計算平均形狀到外輪廓的變換T：

(2b)基于小鼠體外輪廓與腎臟之間的位置相關(guān)性，將配準變換T近似為腎臟平均形狀到待分割圖像腎臟的變換，對腎臟平均形狀做變換T，得到腎臟初始位置P^O，即

(2c)根據(jù)步驟(2b)中確定的腎臟初始位置P^O，使用步驟(2a)所述方法分割得到腎臟初始分割結(jié)果T^O；

(2d)對步驟(2c)得到的腎臟輪廓T^O的點云數(shù)據(jù)使用三角剖分算法進行面片化，再對得到的三角網(wǎng)格數(shù)據(jù)通過Mesh voxelization算法進行體素化。

步驟4：初始化Graph

(3a)Graph的構(gòu)造：

待分割圖像中的每一個像素點抽象為Graph中的頂點，每個頂點p有兩種邊:

③t-link,邊e_{p,S}，e_{p,T}分別表示像素點和前景(腎臟)、背景(非腎臟)的

關(guān)系，權(quán)值為

R(L)+k₁S_p(L)

其中

R(L)＝-ln Pr(I_r|L)

其中P_r(I_r|L)使用步驟(1)中由灰度值得到的分類函數(shù)計算；

當(dāng)時，若點p在內(nèi),否則

當(dāng)時，Pr_p(L)使用步驟(1)由坐標得到的分類函數(shù)計算，T_n為分類函數(shù)決策邊界處像素點坐標與訓(xùn)練樣本分割邊界的平均距離；

④n-link，表示像素點和8鄰居系統(tǒng)中其他八個像素點之間的關(guān)系，權(quán)值為

B(L)+k₂S_p1,p2,

其中

σ為經(jīng)驗值；

其中

L_p1，L_p2為點p₁，p₂的標簽(腎臟或非腎臟)；為點p_i到平均形狀 的歐氏距離；L為訓(xùn)練圖像的數(shù)目，σ_i²由第i個訓(xùn)練圖像計算 得到，N_i為第i個訓(xùn)練圖像中像素點的個數(shù)， 為點p_j到平均形狀的歐氏距離；

(3b)使用步驟(2)得到的預(yù)分割結(jié)果對步驟(3a)構(gòu)造的Graph進行初始化:根據(jù)步驟(2)預(yù)分割結(jié)果，在腎臟內(nèi)部的像素點所對應(yīng)的Graph中頂點的n-link置為e_{p,S}＝Inf,e_{p,T}＝0；在腎臟外部的像素點對應(yīng)的n-link置為e_{p,S}＝0,e_{p,T}＝Inf，得到初始化后的圖G_l。

步驟5：使用Boykov-Kolmogorov的Graph Cut算法對步驟(3)中得到的初始化后的圖G_l計算其最大流，即最小化能量方程

E(L)＝λ(R(L)k₁S_p(L))+(B(L)+k₂S_p1,p2)，

得到分割G_l并使得能量方程取值最小的一個割(Cut),即得到最終的腎臟邊界。

以上所述，僅為本發(fā)明中的一個具體實例，但本發(fā)明的保護范圍并不局限于此，任何熟悉該技術(shù)的人在本發(fā)明所揭露的技術(shù)范圍內(nèi)，可理解想到的變換或替換，都應(yīng)涵蓋在本發(fā)明的包含范圍之內(nèi)，因此，本發(fā)明的保護范圍應(yīng)該以權(quán)利要求書的保護范圍為準。