基于配準(zhǔn)分割的G raph Cuts自動(dòng)分割算法在肝臟圖像中的研究

王建軍　謝勤嵐

（中南民族大學(xué)生物醫(yī)學(xué)工程學(xué)院　武漢　430074）

1　引言

肝癌是最常見的6種癌癥之一，同時(shí)也是全世界范圍內(nèi)導(dǎo)致人類死亡的主要原因。在2012年，全世界大約有78萬人被診斷為肝癌，約有74萬人死于肝癌［1］。因此，較早的診斷和治療肝病是非常重要的。肝臟醫(yī)學(xué)影像分割非常重要，它可以輔助醫(yī)生對(duì)肝臟疾病患者的肝臟進(jìn)行診斷、功能評(píng)定和治療決策。但是在醫(yī)學(xué)圖像處理中，肝臟的分割一直都是一項(xiàng)巨大挑戰(zhàn)的任務(wù)。如圖1所示的肝臟解剖結(jié)構(gòu)，其周圍有主動(dòng)脈、骨頭、腎等灰度相似的器官，還有腫瘤等疾病，肝臟對(duì)比度低、邊界模糊，很難將其完全分割出來。

圖1　肝臟的解剖結(jié)構(gòu)

肝臟分割方法［2］包括半自動(dòng)分割和全自動(dòng)分割兩類。目前已經(jīng)有許多分割算法能對(duì)肝臟進(jìn)行有效的分割。周向榮等使用概率模型估計(jì)初始空間的位置并計(jì)算肝臟分布的概率［3］，能夠?qū)崿F(xiàn)對(duì)肝臟的全自動(dòng)分割，但分割精度不夠理想。Rikxoort等通過使用K最近鄰域方法計(jì)算肝臟中每個(gè)像素的概率和多尺度Atlas配準(zhǔn)方法［4］，對(duì)肝臟進(jìn)行有效的全自動(dòng)分割。Foruzan等提出基于灰度分析和結(jié)構(gòu)信息的方法［5］，使用期望最大值算法，計(jì)算肝臟灰度范圍的統(tǒng)計(jì)參數(shù)，結(jié)合閾值方法，能夠?qū)⒏闻K與周圍的組織有效的分割出來。R Pohle，KD Toennies提出的模型自適應(yīng)區(qū)域增長(zhǎng)算法［6］，對(duì)于對(duì)比度不高的肝臟CT影像，都能獲得一個(gè)較好的分割結(jié)果。但是他們的方法對(duì)于初始種子點(diǎn)的選取非常敏感，魯棒性不夠好。

本文利用基于配準(zhǔn)分割的Graph Cuts自動(dòng)分割算法，對(duì)人體肝臟CT圖像進(jìn)行自動(dòng)分割。該方法的分割流程如圖2所示，主要的思想是利用Rikxoort提出的基于Atlas的分割方法對(duì)肝體進(jìn)行有效的分割，利用權(quán)重投票［7］將所有的Atlas的分割結(jié)果融合成一幅預(yù)分割圖像，然后利用水平集分割中的符號(hào)距離函數(shù)［8］自動(dòng)標(biāo)記預(yù)分割圖像目標(biāo)和背景的像素，最后將標(biāo)記的像素標(biāo)簽送入Graph Cuts圖割算法［9］中，完成自動(dòng)分割。

圖2　自動(dòng)分割流程圖

2　基于A tlas的分割

Atlas分割算法的主要思想是利用已經(jīng)分割的Atlas圖像對(duì)待分割圖像進(jìn)行非剛性配準(zhǔn)，通過配準(zhǔn)得到的形變場(chǎng)參數(shù)應(yīng)用到Atlas的二值圖像上，該二值圖像通過應(yīng)用形變場(chǎng)參數(shù)得到的圖像即為待分割圖像的分割結(jié)果。也就是說，Atlas分割算法是將分割問題轉(zhuǎn)化為了配準(zhǔn)的問題。圖3是基于Atlas的分割框架圖。

圖3　Atlas分割框架

Atlas分割算法對(duì)肝臟圖像分割結(jié)果的優(yōu)劣，主要依靠配準(zhǔn)的質(zhì)量。對(duì)于配準(zhǔn)，主要分為四個(gè)步驟：配準(zhǔn)框架（Registration framework）、度量方法（Metrics），插值（Interpolators）和優(yōu)化（Optimisers）。

本文的配準(zhǔn)框架是使用剛性的全局仿射粗配和非剛性的局部B樣條精配組合的框架，其框架模型的公式如式（1）。Igabol(x，y，z) 是配準(zhǔn)的全局仿射粗配模型，Ilocal(x，y，z)是局部精配B樣條自由形變模型（Free-form deformation，F(xiàn)FD）。

度量方法是配準(zhǔn)中較為重要的部分，它是衡量?jī)煞鶊D像配準(zhǔn)優(yōu)劣的標(biāo)準(zhǔn)。目前度量的方法有：均方差（Mean squared difference，MSD），歸一化相關(guān)系數(shù)（Normalised correlation coefficient，NCC）［10］，歸一化互信息（Normalized Mutual Information，NMI），Kappa統(tǒng)計(jì)（Kappa statistic，KS）［11］以及本文所采用的互信息（Mutual information，MI）［12］。MI的公式定義如式（2）。

其中，IF是固定圖像，IM是浮動(dòng)圖像，μ是形變參數(shù)，LF、LM分別是固定圖像和浮動(dòng)圖像的灰度級(jí)，pF和 pM是邊緣概率，p( f ，m;μ )是兩幅圖像的聯(lián)合概率密度。

圖像插值的作用主要是在經(jīng)過變換的像素位置可能不再是整數(shù)后獲得整數(shù)位置的像素值，優(yōu)化是是變換參數(shù)的相似性測(cè)度達(dá)到最優(yōu)。

3　權(quán)重投票

在Atlas分割實(shí)驗(yàn)中，對(duì)于待分割圖像有多個(gè)Atlas圖像，就有多少個(gè)分割結(jié)果。然而，這些分割結(jié)果的精度差別較大，并不是所有的分割結(jié)果都適合做Graph Cuts實(shí)驗(yàn)的預(yù)分割。為了得到適合Graph Cuts分割實(shí)驗(yàn)的預(yù)分割圖像，可以采用聯(lián)合所有的Atlas圖像形變分割融合成一幅分割結(jié)果圖像。權(quán)重投票融合公式［13］如下其中，c是分類數(shù)，本算法只有兩類即目標(biāo)和背景，故 c為2；lc(x)是Atlas分割結(jié)果中，圖像的像素x屬于目標(biāo)或者背景的概率；δ[?]是克羅內(nèi)克函數(shù)，當(dāng)c和(Li?Ti)(x)相等時(shí)，其值為1，當(dāng) c和(Li?Ti)(x)不相等時(shí)，其值為0；wi是標(biāo)量權(quán)重因子，它是在Atlas試驗(yàn)中，從配準(zhǔn)的文件中提取的NMI的值。

權(quán)重投票融合，可以將多個(gè)Atlas分割結(jié)果融合為一幅圖像，將融合后的圖像作為Graph Cuts分割實(shí)驗(yàn)的預(yù)分割圖像，可以增加實(shí)驗(yàn)的魯棒性。

4　水平集符號(hào)距離函數(shù)

水平集函數(shù)是一種跟蹤輪廓和曲面演化的數(shù)值方法，而不是直接計(jì)算曲面的輪廓。如圖4所示，是零水平集函數(shù)示意圖。圖中黑色粗線輪廓部分，是零水平集函數(shù)代表的輪廓，它是由更高維度的水平集函數(shù)的等于0時(shí)對(duì)應(yīng)的曲線輪廓。

圖4　零水平集函數(shù)示意圖

水平集符號(hào)距離函數(shù)表示像素點(diǎn)到零水平集函數(shù)的歐式距離。在本文的實(shí)驗(yàn)中，主要是像素點(diǎn)到Atlas分割結(jié)果圖像中目標(biāo)的邊緣的距離，距離邊緣越遠(yuǎn)的像素，距離越大，距離邊緣越近的像素，距離越小。本文符號(hào)距離函數(shù)的完成是基于作者C.R.Maurer提出的線性時(shí)間算法［14］。圖5是利用線性時(shí)間算法求Atlas分割結(jié)果圖像的符號(hào)距離函數(shù)圖像的過程，（a）是Atlas分割完成的二值圖像，（b）是線性時(shí)間算法計(jì)算的符號(hào)距離函數(shù)圖像，從圖中可以看出，距離二值圖像邊緣越遠(yuǎn)的像素點(diǎn)灰度值越亮，表明該到邊緣的歐式距離越遠(yuǎn)。

圖5　Atlas分割結(jié)果和符號(hào)距離函數(shù)圖像

圖6是利用符號(hào)距離函數(shù)標(biāo)記目標(biāo)和背景像素的過程。圖（a）是權(quán)重投票的結(jié)果，將所有的Atlas分割結(jié)果融合為一幅預(yù)分割圖像，圖（b）是使用三種不同大小的符號(hào)距離標(biāo)記融合圖像的目標(biāo)和背景，其中的內(nèi)符號(hào)距離3的大小決定著白色目標(biāo)像素的大小，外符號(hào)距離1和外符號(hào)距離2決定著背景像素大小。

圖6　利用符號(hào)距離函數(shù)標(biāo)記前背景像素

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

本文所有的實(shí)驗(yàn)都是基于IKT（InsightSegmentation and Registration Toolkit）［15］框架醫(yī)學(xué)圖像處理工具。其中Atlas分割模塊是采用開源配準(zhǔn)軟件包Elastix，權(quán)重投票和水平集符號(hào)距離函數(shù)等模塊是采用ITK的內(nèi)部類，Graph Cuts算法的基礎(chǔ)代碼部分是由 Tustison N.和 Yushkevich P.［16］等提供。本實(shí)驗(yàn)的數(shù)據(jù)是來自臨床的CT圖像數(shù)據(jù)集3Dircadb。3Dircadb數(shù)據(jù)集是由20組CT醫(yī)學(xué)影像組成，圖像切片均為512×512像素，切片的像素間隔為0.56~0.87，切片層數(shù)范圍是91~260。

在Atlas分割實(shí)驗(yàn)中，用3Dircadb數(shù)據(jù)集其中一個(gè)病人圖像作為訓(xùn)練集，另外19個(gè)病人圖像作為測(cè)試集，其結(jié)果如表1所示。Patient_ID表示病人CT圖像的序號(hào)，Overlap是Atlas分割結(jié)果圖像融合后的分割精度。圖7是對(duì)3Dircadb數(shù)據(jù)中第一個(gè)病人圖像使用Atlas分割算法后將其利用權(quán)重融合的示意圖。

表1　融合Atlas分割重疊率結(jié)果

圖7　Atlas分割結(jié)果和權(quán)重投票融合的例子

將Atlas分割結(jié)果融合后，需要對(duì)其進(jìn)行像素標(biāo)記。圖8是利用符號(hào)距離函數(shù)標(biāo)記目標(biāo)和背景像素標(biāo)簽的過程，圖8（a）是待分割圖像，圖8（b）是融合Atlas分割結(jié)果圖像，圖8（c）利用符號(hào)距離函數(shù)標(biāo)記像素的前背景像素圖像。在圖8（c）中，藍(lán)色部分是融合Atlas分割結(jié)果，白色部分是標(biāo)記目標(biāo)像素，藍(lán)色外的灰白色部分是標(biāo)記的背景像素。經(jīng)過大量的實(shí)驗(yàn)發(fā)現(xiàn)，基于Atlas分割結(jié)果圖像，當(dāng)內(nèi)符號(hào)距離為大于5的像素全部標(biāo)記為目標(biāo)，外符號(hào)距離為0～15的像素標(biāo)記為背景時(shí)Graph Cuts自動(dòng)分割效果最好，其結(jié)果要優(yōu)于融合Atlas分割結(jié)果圖像的分割精度。

圖8　利用符號(hào)距離函數(shù)標(biāo)記像素

圖9是對(duì)數(shù)據(jù)集3Dircadb中某一個(gè)病人的采用Atlas分割算法和Graph Cuts自動(dòng)分割算法的對(duì)比圖，圖10是兩種算法對(duì)數(shù)據(jù)集20個(gè)病人分割結(jié)果的線箱圖。從圖中可以看出，Graph Cuts自動(dòng)分割算法不盡在割精度上優(yōu)于Atlas分割算法，而且魯棒更好。

圖9　兩種分割結(jié)果對(duì)比圖

圖10　兩種算法分割結(jié)果重疊率Box-plot圖

6　結(jié)語

本文針對(duì)肝臟CT醫(yī)學(xué)影像的分割，采用Atlas配準(zhǔn)分割為預(yù)分割結(jié)果，通過權(quán)重投票融合預(yù)分割結(jié)果圖像，最后使用水平集符號(hào)距離函數(shù)標(biāo)記前背景像素，實(shí)現(xiàn)GraphCut自動(dòng)分割。最終的實(shí)驗(yàn)結(jié)果表明，本文提出的基于配準(zhǔn)分割的Graph Cuts自動(dòng)分割算法，相較于傳統(tǒng)的Atlas配準(zhǔn)分割有較大精度的提升，在操作性方面，優(yōu)于傳統(tǒng)的Graph Cuts分割算法，實(shí)現(xiàn)自動(dòng)分割目標(biāo)。

［1］Ferlay J，Soerjomataram I，Dikshit R，et al.Cancer incidence andmortality worldwide：sources，methods andmajor patterns in GLOBOCAN 2012［J］.International JournalofCancer，2015，136（5）：359-363.

［2］Mharib A M，Ram li A R，Mashohor S，et al.Survey on liver CT image segmentationmethods［J］.Artificial Intelligence Review，2012，37（2）：83-95.

［3］Zhou X，Kitagawa T，Hara T，etal.Constructing a Probabilistic Model for Automated Liver Region Segmentation Using Non-contrast X-Ray Torso CT images［J］.Med Image ComputComputAssist Interv，2006，9（2）：856-863.

［4］Rikxoort EMV，Arzhaeva Y，Ginneken BV.Automatic segmentation of the liver in computed tomograpy scans with voxel classification and atlasmatching［J］.3d Segmentation in the Clinic A Grand Challenge，2007，147（6）：132-136.

［5］Foruzan A H，Zoroofi R A，HoriM，et al.Liver segmentation by intensity analysis and anatomical information in multi-slice CT images［J］.International Journal of Computer Assisted Radiology and Surgery，2009，4（3）：287-297.

［6］ Pohle R，Toennies K D.A New Approach for Model-Based Adaptive Region Growing in Medical Image Analysis［M］.Computer Analysis of Images and Patterns.Springer Berlin Heidelberg，2001：238-246.

［7］Rohlfing T，Jr C R M.Multi-classifier framework for atlas-based image segmentation［J］.Pattern Recognition Letters，2005，26（13）：2070-2079.

［8］J.A.Sethian.Level SetMethods and FastMarching Methods［M］.Cambridge University，Second edition，1999，7（2）：81-85.

［9］Boykov Y Y，Jolly MP.Interactive Graph Cuts for Optimal Boundary&Region Segmentation of Objects in N-D Images［J］.Iccv，2001，1：105-112.

［10］Knops Z F，Maintz JB A，Viergever MA，et al.Normalized Mutual Information Based PET-MR Registration Using K-Means Clustering and Shading Correction［J］.Medical Image Analysis，2006，10（3）：432-439.

［11］Fitzpatrick JM，Hill D LG，Shyr Y，etal.Visualassessment of the accuracy of retrospective registration of MR and CT images of the brain［J］.IEEE Transactions on Medical Imaging，1998，17（4）：571-85.

［12］Stefan Klein，Marius Staring.Elastix themanual［M］.Image Sciences Institute，2015，9（4）：5-6.

［13］Klein S，Van-Der-Heide U，Lips I，et al.Automatic segmentation of the prostate in 3D MR images by atlas matching using localizedmutual information［J］.Medical Physics，2008，35（4）：1407-1417.

［14］Jr CR M，QiR，Raghavan V.A Linear Time Algorithm for Computing Exact Euclidean Distance Transforms of Binary Images in Arbitrary Dimensions［J］.Pattern Analysis&Machine Intelligence IEEE Transactions on，2003，25（2）：265-270.

［15］Ibanez L，Schroeder W，Ng L，et al.The ITK Software Guide［J］.Computational Statistics&Data Analysis，2005，21：231-256.

［16］N.J TUSTISON，P.Yushkevich，Z.Song.Graph Cuts，Caveat Utilitor，and Euler's Bridges of Konigsberg［J］.Insight Journal，2008，11（14）：128-138.