趙潔，蔣世忠，黃展鵬，歐陜興

（1.廣東藥科大學(xué)醫(yī)藥信息工程學(xué)院，廣州510006；2.南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院廣東省醫(yī)學(xué)圖像處理重點(diǎn)實(shí)驗(yàn)室，廣州510515；3.廣州軍區(qū)總醫(yī)院放射科，廣州510010）

1 引 言

據(jù)世界衛(wèi)生組織2008年公布的統(tǒng)計(jì)結(jié)果，冠心病是當(dāng)前世界范圍內(nèi)引發(fā)死亡的主要原因[1]。早診斷、早治療是降低死亡率的方法。雙源CT圖像采集速度比多層螺旋CT快，可以在一個(gè)心動(dòng)周期內(nèi)采集心臟圖像，并且輻射劑量要減少一半以上，在冠心病診斷中發(fā)揮著重要作用。但是，一名患者的心臟CT一般要產(chǎn)生大于400張圖像，讀片會(huì)耗費(fèi)醫(yī)生大量精力，增加誤判率。如何幫助醫(yī)生從大量CTA數(shù)據(jù)中準(zhǔn)確快速的分割出冠狀動(dòng)脈，是我們研究的主要問題[2]。

冠狀動(dòng)脈的快速準(zhǔn)確分割是血管定量、定性測量和疾病診斷的前提。但是由于心臟CTA圖像數(shù)據(jù)量大、冠狀動(dòng)脈和其周邊組織灰度比較相近、運(yùn)動(dòng)偽影和外部磁場引起的噪聲干擾等問題，準(zhǔn)確分割非常困難。目前，國內(nèi)外對冠狀動(dòng)脈提取常規(guī)的方法分有：區(qū)域生長法、統(tǒng)計(jì)學(xué)方法、活動(dòng)輪廓模型法、中心線法、多尺度模型，下面是幾種方法的介紹[3-4]。

區(qū)域生長法較為簡單，速度快，但是生長條件和種子點(diǎn)需人為指定，由于冠狀動(dòng)脈與心室灰度相近，當(dāng)受到噪聲和對比劑不均勻影響時(shí)，目標(biāo)局部會(huì)出現(xiàn)孔洞和泄露問題。Metz等[5]在經(jīng)典區(qū)域生長算法的基礎(chǔ)上加入分支和泄露檢測，提取出冠狀動(dòng)脈的三個(gè)大的血管分支，但是小分支仍無法分割。楊榮騫等[6，7]用形態(tài)學(xué)方法和區(qū)域生長法從雙源CT數(shù)據(jù)中分割出心腔和冠狀動(dòng)脈。Tuchschmid等[8]用波傳播技術(shù)，通過與空間結(jié)構(gòu)一致的方向生長出冠狀動(dòng)脈。

統(tǒng)計(jì)學(xué)方法包括聚類法、分類器法、馬爾科夫模型等。Yan Yang在2004年IEEE EMBS國際會(huì)議[9]提出一種基于CT的冠狀動(dòng)脈三維分割算法，具有一定代表性。Lacoste C等[10]利用Markov標(biāo)記點(diǎn)進(jìn)行2D造影圖像的冠狀動(dòng)脈分割，但是，目前該方法僅限于2D數(shù)據(jù)，由于3D造影數(shù)據(jù)計(jì)算量龐大，尚未能擴(kuò)展。

活動(dòng)輪廓法在醫(yī)學(xué)圖像分割上用的較多，效果較好。1988年，Osher和 Sethian[11]提出了利用水平集方法來求解幾何曲線的演化。水平集的演化曲線是閉合的，可以較好的收斂到目標(biāo)邊緣。但是，傳統(tǒng)水平集存在對噪聲敏感，弱邊緣不易檢測，易越過邊緣的問題[12]。在對冠狀動(dòng)脈分割上，選取的初始輪廓必須在目標(biāo)附近，并且由于和鄰近組織灰度相近，容易擴(kuò)展到心腔。針對血管泄露問題，Yan Yang等[13]通過在水平集中加入統(tǒng)計(jì)信息進(jìn)行了改進(jìn)，可以分割出左冠狀動(dòng)脈，左前降支，左旋支和右冠狀動(dòng)脈的主干，但是更細(xì)的血管未實(shí)現(xiàn)分割。

中心線血管分割主要包括：（1）跟蹤法，（2）最小路徑技術(shù)。跟蹤法通常需要在初始點(diǎn)和分支點(diǎn)進(jìn)行人工交互，為了自動(dòng)判斷分支結(jié)構(gòu)，Piasque[14]提出利用球和平行六面體的連通分量來判斷，蔣世忠、趙潔[15]和黃展鵬[16]提出構(gòu)造球模型，通過層次聚類算法來自動(dòng)判斷分支。最小路徑技術(shù)有利于局部病變血管的提取，快速行進(jìn)法是最小路徑技術(shù)的優(yōu)化，可通過優(yōu)化計(jì)算代價(jià)函數(shù)的歐幾里得距離，產(chǎn)生亞像素精度的血管路徑[5]。

通過在不同尺度下分析Hessian矩陣的特征值，可以判斷每一個(gè)體素點(diǎn)是否屬于管狀結(jié)構(gòu)。基于多尺度空間，F(xiàn)rangi[17]提出血管增強(qiáng)方法，Manniesing[18]提出血管增強(qiáng)擴(kuò)散模型，均可以有效增強(qiáng)圖像中的管狀結(jié)構(gòu)。

由于心臟CTA圖像數(shù)據(jù)量大、冠狀動(dòng)脈和其周邊組織灰度比較相近、運(yùn)動(dòng)偽影和外部磁場會(huì)引起噪聲干擾等因素，血管的準(zhǔn)確分割非常困難，而錯(cuò)誤的分割會(huì)使重建和測量精確度受到影響，影響醫(yī)生診斷。在已有的心臟血管分割文獻(xiàn)中，一些方法在二維圖像中可以很好的分割出血管，但是模型復(fù)雜、速度較慢，并不適合三維分割。本研究提出的三維自動(dòng)分割方法，考慮了血管的形態(tài)、灰度和鄰域特征，基于Hessian矩陣、統(tǒng)計(jì)模型和26鄰域生長算法，可以自動(dòng)分割冠狀動(dòng)脈并精確測量血管狹窄信息。

2 算法設(shè)計(jì)

為了檢測血管，本算法同時(shí)從血管形態(tài)、灰度特征和鄰域關(guān)系三個(gè)方面考慮。血管的形態(tài)是管狀的，通過對每個(gè)像素點(diǎn)求Hessian矩陣的特征值增強(qiáng)管狀結(jié)構(gòu)。增強(qiáng)后的圖像，血管區(qū)域灰度偏亮，背景區(qū)域偏暗，利用統(tǒng)計(jì)模型，對血管類和背景類進(jìn)行高斯建模。在血管類中，通過26鄰域生長來約束鄰域關(guān)系，生長出冠狀動(dòng)脈。算法的流程見圖1。

圖1 算法流程圖Fig 1 the flow chart of the proposed method

3 Hessian矩陣

Hessian矩陣可用來區(qū)分圖像中的管狀，片狀，點(diǎn)狀結(jié)構(gòu)。冠狀動(dòng)脈的形態(tài)是近似管狀，因此，利用Hessian矩陣的多尺度濾波可去除其他結(jié)構(gòu)，增強(qiáng)血管。首先，對每個(gè)體素做Hessian矩陣運(yùn)算，見公式（1）。

Hessian矩陣的差分被定義為圖像與δ階高斯倒數(shù)的卷積，見公式（2）。

對每個(gè)體素點(diǎn)的Hessian矩陣計(jì)算特征值λ1，λ2和λ3。特征值判斷準(zhǔn)則見公式（3），可以區(qū)分管狀結(jié)構(gòu)，如血管[15]。

Frangi定義了血管檢測的方程：

其中：

為增強(qiáng)冠狀動(dòng)脈，消除其他組織和主動(dòng)脈，濾波尺度定為[13]，步長為0.5。

4 概率決策分類

高斯混合模型假設(shè)每一類組織符合同一高斯分布，整個(gè)體數(shù)據(jù)的分布由各組織類的高斯概率線性組合而成。

在上一步驟中，通過Hessian矩陣增強(qiáng)了圖像中的管狀結(jié)構(gòu)。假設(shè)增強(qiáng)后的圖像分為2類，分別為背景和管狀區(qū)。管狀區(qū)包含血管、骨骼等管狀結(jié)構(gòu)和少量噪聲?？偟捏w數(shù)據(jù)概率分布由這2類組織的高斯概率密度混合而成，見式（5），wk為每一類所占的權(quán)重。

每一類組織符合同一高斯分布，見式（6）。

參數(shù)wk，uk，σk依靠最大期望算法估計(jì)。

E步利用分布參數(shù)計(jì)算每個(gè)實(shí)例的聚類概率，見式（7），t代表迭代步數(shù)。

M步用最大似然估計(jì)重新估計(jì)分布參數(shù)，得到：

EM兩步反復(fù)迭代直至收斂，得到高斯混合模型的參數(shù)。然后對每一個(gè)體素判別屬于哪一類，判別條件為 wkp（x｜k）＞wjp（x｜j），k～＝j(luò)時(shí)屬于第 k類[19-20]。

5 三維區(qū)域生長

概率決策分類后數(shù)據(jù)被分為2類，血管類和背景[21]。在冠狀動(dòng)脈上設(shè)置種子點(diǎn)，在三維體數(shù)據(jù)中進(jìn)行區(qū)域生長，生長步驟如下：

（1）將種子點(diǎn)標(biāo)記為1。

（2）搜索種子點(diǎn)26鄰域，如果閾值小于設(shè)定值，則合并到種子區(qū)域，并將其置為1。

（3）直到無鄰域點(diǎn)滿足條件時(shí)生長截止。

由于在三維區(qū)域生長算法之前約束了血管形狀和灰度，去除了心室等區(qū)域，可以在生長過程中避免泄漏。

6 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)數(shù)據(jù)來自廣州軍區(qū)總醫(yī)院的SOMATOM Definition雙源CT，采用心電門控技術(shù)獲取的心臟雙源CT三維數(shù)據(jù)，分辨率為512×512×411。圖2為位于一組CT序列的第1幀、第40幀和第60幀切片。

圖2 心臟CTA圖像（第1、40、60幀）Fig 2 Cardiac Computed tomography angiography images（slice1，slice40，slice60）

從圖2第40幀和第60幀切片可以看出，左前降支與左回旋支均和左心室空間相鄰并且灰度接近。如果直接使用區(qū)域生長或活動(dòng)輪廓等方法進(jìn)行分割，會(huì)發(fā)生灰度泄露。

用本研究提出的新方法，首先對體數(shù)據(jù)每個(gè)像素進(jìn)行Hessian運(yùn)算，利用式（4）檢測出血管，得到圖3，可以看出管狀結(jié)構(gòu)被增強(qiáng)。

圖3 多尺度濾波結(jié)果（第1、40、60幀）Fig 3 Multi-scale filtering results（slice1，slice40，slice60）

用雙高斯模型，多尺度濾波后的數(shù)據(jù)可分為背景和血管兩大類，因此，三維的血管分割問題變?yōu)榱穗p高斯模型的參數(shù)估計(jì)問題。根據(jù)EM算法，估計(jì)出每一類組織的均值和方差，擬合出圖4。

圖4 雙高斯直方圖擬合Fig 4 Double Gaussian histogram fitting

以圖4中的參數(shù)進(jìn)行分割，得到分類結(jié)果，黑色為背景，白色為血管，見圖5。

圖5 統(tǒng)計(jì)分類結(jié)果（第1、40、60幀）Fig 5 Statistical classification results（slice1，slice40，slice60）

在冠狀動(dòng)脈與升主動(dòng)脈結(jié)合部分任選1點(diǎn)作為種子點(diǎn)，進(jìn)行26鄰域區(qū)域生長，生長結(jié)果見圖6。因?yàn)榉诸惤Y(jié)果已被處理為黑白圖，所以生長速度非?？?。

圖6 三維區(qū)域生長結(jié)果（第1、40、60幀）Fig 6 Three dimensional region growing results（slice1，slice40，slice60）

將分割好冠狀動(dòng)脈進(jìn)行三維重建，可以看到血管的三維效果圖，見圖7，左冠狀動(dòng)脈重建效果較好，避免了血管泄露問題，無偽血管，動(dòng)脈狹窄清晰可見。

圖7 左冠狀動(dòng)脈三維重建結(jié)果Fig 7 Three dimensional reconstruction of left coronary artery

圖8 中心線提取結(jié)果Fig 8 Centerline extraction result

圖8中明顯有冠狀動(dòng)脈狹窄，狹窄處僅為1～2個(gè)像素寬，但本文算法依然可以正確分割，沒有發(fā)生斷裂。

為了對所提算法進(jìn)行定量分析，本研究測試了MICCAI08中已手工標(biāo)注了冠狀動(dòng)脈中心線的數(shù)據(jù)，數(shù)據(jù)來源于SOMATOM Definition，分辨率為512×512×288。圖9為算法分割的冠狀動(dòng)脈中心線與金標(biāo)準(zhǔn)之間的對比。

根據(jù)鹿特丹冠狀動(dòng)脈評(píng)價(jià)軟件計(jì)算[22]，本研究分割的中心線與金標(biāo)準(zhǔn)之間的重合率為89%，經(jīng)廣州軍區(qū)總醫(yī)院放射科醫(yī)生評(píng)價(jià)，主要的三級(jí)冠狀動(dòng)脈基本提取，符合醫(yī)學(xué)診斷要求。

7 小結(jié)

本研究所提算法結(jié)合了血管灰度特征、形態(tài)特征和結(jié)構(gòu)特征，基于統(tǒng)計(jì)分類，多尺度Hessian血管增強(qiáng)和三維區(qū)域生長構(gòu)建出心臟血管模型。新算法對冠狀動(dòng)脈的分割比較準(zhǔn)確，避免了灰度泄露等問題，重建效果好，能很好的應(yīng)用于心臟雙源CTA序列圖像。正確的分割出冠狀動(dòng)脈是心血管疾病診斷的基礎(chǔ)，下一步工作要解決的問題是自動(dòng)檢測冠狀動(dòng)脈狹窄，輔助醫(yī)生實(shí)施合理的心血管疾病輔助診斷和心臟手術(shù)規(guī)劃。

圖9 本文算法中心線分割結(jié)果與金標(biāo)準(zhǔn)比較（a）本文算法；（b）金標(biāo)準(zhǔn)Fig 9 Comparison between the results of the proposed method and golden standard（a）the proposed algorithm（b）golden standard