周應(yīng)杰，吳海堅(jiān)

廣州醫(yī)學(xué)院附屬腫瘤醫(yī)院， 廣東 廣州 510095

基于偏微分方程的CT圖像分割及其MATLAB實(shí)現(xiàn)

周應(yīng)杰，吳海堅(jiān)

廣州醫(yī)學(xué)院附屬腫瘤醫(yī)院， 廣東 廣州 510095

CT圖像在醫(yī)學(xué)影像領(lǐng)域的應(yīng)用最為廣泛，CT圖像的分割對(duì)病變區(qū)域的提取、組織定位、組織測量以及實(shí)現(xiàn)三維重建有著非常重要的作用。近年來，基于偏微分方程的圖像分割作為一種比較新穎且有效的圖像分割方法，逐漸成為研究的熱點(diǎn)。本文對(duì)基于偏微分方程(PDEs)的GAC和C-V兩種圖像分割模型進(jìn)行了探討，利用MATLAB語言編程實(shí)現(xiàn)，并結(jié)合實(shí)際CT圖像進(jìn)行了實(shí)驗(yàn)研究，為這些模型在CT圖像的分割上的應(yīng)用提供科學(xué)依據(jù)。

CT圖像；CT圖像分割；偏微分方程；GAC模型；C-V模型

0 前言

圖像分割一直是圖像處理研究中的重要任務(wù)和熱點(diǎn)問題，它早已被廣泛地應(yīng)用于醫(yī)學(xué)圖像處理之中[1]。近十幾年來，隨著偏微分方程方法在圖像處理中應(yīng)用的興起，圖像分割的偏微分方程方法也逐步取得了很大的進(jìn)展。圖像分割可以理解為將圖像中有意義的特征區(qū)域或者需要應(yīng)用的特征區(qū)域提取出來，這些特征區(qū)域可以是像素的灰度值、物體輪廓特性曲線、紋理特性等，也可以是空間頻譜或直方圖特征等[2]。醫(yī)學(xué)圖像分割是醫(yī)學(xué)圖像處理和分析中的關(guān)鍵技術(shù)，分割后的圖像被廣泛應(yīng)用于組織的定量分析、診斷、病變組織的定位、解剖結(jié)構(gòu)的學(xué)習(xí)、治療規(guī)劃、功能成像數(shù)據(jù)的局部體效應(yīng)校正和計(jì)算機(jī)指導(dǎo)手術(shù)等[3]。隨著影像醫(yī)學(xué)在臨床醫(yī)學(xué)的成功應(yīng)用，圖像分割在影像醫(yī)學(xué)中發(fā)揮著越來越大的作用[4]。CT圖像是利用計(jì)算機(jī)技術(shù)處理人體組織器官得到的切面顯像，它可清晰地顯示人體內(nèi)的組織結(jié)構(gòu)。CT技術(shù)作為一種先進(jìn)的疾病診斷手段，廣泛應(yīng)用于醫(yī)學(xué)領(lǐng)域中，但是由于醫(yī)學(xué)圖像具的復(fù)雜性和多樣性，使得對(duì)CT圖的分割很困難。目前很多學(xué)者致力于此方面的研究，取得了很多成果，但至今還沒有找到一種十分滿意的方法應(yīng)用于臨床醫(yī)療[5]。

1 傳統(tǒng)的圖像分割方法和基于PDE的圖像分割模型

在傳統(tǒng)的圖像處理技術(shù)中，已經(jīng)有相當(dāng)多的研究成果和方法，這些方法可以歸為3大類：基于區(qū)域的分割方法、基于邊緣的分割方法和基于閾值的分割方法[6-7]?；陂撝档姆指罘椒ㄏ鄬?duì)簡單，計(jì)算量小，穩(wěn)定性較好，但抗干擾能力不強(qiáng)，對(duì)目標(biāo)和背景灰度有梯度變化的圖象效果較差?；趨^(qū)域的分割方法，依賴于圖像的空間局部特征，如灰度、紋理及其他像素統(tǒng)計(jì)特性的均勻性等。典型的基于區(qū)域的分割方法有區(qū)域生長、區(qū)域分裂以及區(qū)域生長與分裂相結(jié)合的方法等，由于這些方法直接依賴于圖像的灰度值，因此它們的主要優(yōu)點(diǎn)是對(duì)噪聲不敏感，但這些方法常會(huì)造成圖像的過分分割問題，并且分割結(jié)果較大程度上依賴于種子點(diǎn)的選擇，分割所得到的區(qū)域的形狀也依賴于所選擇的分割算法?；谶吔绲姆指罘椒ㄖ饕抢锰荻刃畔⒋_定目標(biāo)的邊界，包括局部微分算子、Roberts算子、Prewitt梯度算子、canny算子和Laplacian二階差分算子等。這些方法不依賴于已處理像素的結(jié)果，但缺點(diǎn)是對(duì)噪聲敏感，而且當(dāng)邊緣像素值變化不明顯時(shí)，容易產(chǎn)生假邊界或不連續(xù)的邊界[8]。針對(duì)不同領(lǐng)域的應(yīng)用，相繼又提出了基于統(tǒng)計(jì)學(xué)的方法、基于小波變換的方法、基于信息論的方法、基于模糊分割的方法、基于神經(jīng)網(wǎng)絡(luò)的方法和基于遺傳算法的方法等，每種方法都有它的局限性。

為了克服傳統(tǒng)的圖象處理方法的缺陷，近年來人們提出了基于偏微分方程的圖象分割方法，此方法主要是采用活動(dòng)輪廓模型來實(shí)現(xiàn)圖像分割[9-10]。1987年Kass M等人共同發(fā)表了題為“Snake：Active Contours”的論文[11]，首次提出了運(yùn)用活動(dòng)輪廓模型進(jìn)行2D和3D圖像分割的思想。該模型得到迅速發(fā)展，逐漸成為最有活力和最成功的圖像分割技術(shù)之一。

1.1 測地線活動(dòng)輪廓模型（GAC模型）

活動(dòng)輪廓模型（Active Contours）的基本思想是將圖象分割問題歸結(jié)為最小化一個(gè)封閉曲線C(p)的能量泛函[12]：

根據(jù)曲率運(yùn)動(dòng)演化的性質(zhì)，可以略去第二項(xiàng)，為避免第三項(xiàng)的負(fù)號(hào)，可引進(jìn)一個(gè)新函數(shù)g(r)，r∈R+，(1)式可改寫為：

此模型存在嚴(yán)重缺陷，它不僅依賴于曲線C的幾何形狀和位置，而且還依賴于曲線的參數(shù)p，為了克服此缺陷，Caselles V等人[10]于1997年提出了不含自由參數(shù)的測地線活動(dòng)輪廓（Geodesic active contours，GAC）模型，該模型的提出是PDE（偏微分方程）方法在圖像分割中的重大突破。該模型以最小化以下能量泛函來確定活動(dòng)輪廓：

上式中是閉合曲線C的加權(quán)弧長，由于以上泛函是建立在曲線固有參數(shù)弧長之上的，因而它消除了經(jīng)典蛇模型依賴自由參數(shù)的缺陷?？梢宰C明最小化(3)式對(duì)應(yīng)的梯度下降流為：

即(4)式為GAC模型式。曲線按照此式演化時(shí)將受兩種“力”的支配：① 來自曲線自身曲率運(yùn)動(dòng)的“內(nèi)力”，它受圖像I(x,y)梯度標(biāo)量場g(x,y)的控制，在圖像邊緣這種力變得很小。② 來自于g(x,y)的梯度▽g，由于g(x,y)=g(|▽I(x,y)|)，▽g是由圖像I(x,y)產(chǎn)生的，所以第二項(xiàng)力稱為外力。它能使C向?qū)ο蟮倪吘壙拷?，并穩(wěn)定在邊緣上。從以上分析可知GAC模型存在一個(gè)局限性，即當(dāng)圖像中存在較深的凹陷邊界時(shí)，GAC模型將使C停止在某一能量局部極小值狀態(tài)，并不與對(duì)象邊界一致。此時(shí)若在(4)式中再增加一個(gè)收縮力，使C遇到較深的凹陷時(shí)繼續(xù)收縮，于是GAC模型改造為：

式中c為一可選常數(shù)，(5)式為推廣的GAC模型。

1.2 無邊緣活動(dòng)輪廓模型（C_V模型）

如果圖像既沒有明顯的邊緣，又缺乏明顯的紋理特征，利用以上討論的GAC模型將難以成功的分割。Chan和Vese提出的基于簡化Mumford-Shah模型和水平集方法相結(jié)合[12]可以有效地用于此類圖像分割。其思路是，如果能找到閉合曲線C將圖像分為內(nèi)部Ω1和外部Ω2，使Ω1和Ω2的平均灰度值恰好反映出對(duì)象與背景之間的灰度平均值的差別，那么這一閉合曲線C就可以看作是對(duì)象的輪廓，其能量泛函如下：

上式稱為無邊緣活動(dòng)輪廓模型（Active contours without edges， 或稱C-V模型），其中第一項(xiàng)為C的弧長，第二和第三項(xiàng)為內(nèi)部和外部的灰度值與標(biāo)量C1和C2的平方誤差，只有當(dāng)C達(dá)到正確位置時(shí)這兩項(xiàng)的值才能同時(shí)達(dá)到最小。

采用水平變分集方法，先在(6)式引入Heaviside函數(shù)，將它修改為關(guān)于嵌入函數(shù)u的泛函：

這樣在函數(shù)u固定的條件下相對(duì)C1、C2最小化式(7)，可得：

C1、C2分別是輸入圖像I（x,y）在Ω1（曲線內(nèi)部）和Ω2（曲線外部）的平均值。在C1、C2固定的條件下，相對(duì)于u最小化式(7)可得：

于是通過連立方程(8)和(9)求穩(wěn)態(tài)解，便得到分割結(jié)果。

1.3 基于PDE方法對(duì)CT圖像分割的實(shí)現(xiàn)

1.3.1 GAC模型的數(shù)值解法及實(shí)驗(yàn)結(jié)果

本文采用迎風(fēng)方案實(shí)現(xiàn)推廣GAC模型的水平集方法，限于篇幅，在此不詳細(xì)列出，這里只給出結(jié)果，具體理論方法可以參考微幾何及參考文獻(xiàn)[13-14]。圖1和圖3分別是原始的腦部和體部CT圖，圖2和圖4分別是對(duì)應(yīng)分割后的腦部和體部CT圖，明顯看出，由于GAC模型的特點(diǎn)，此方法最多只能把整個(gè)圖像從背景中分割出來，且分割邊緣也不太理想。

圖1 腦部原始CT圖

圖2 分割后的腦部CT圖

圖3 體部原始圖

圖4 分割后的體部圖

1.3.2 C_V模型的數(shù)值解法及實(shí)驗(yàn)結(jié)果

本文采用半隱式方案實(shí)現(xiàn)變分水平集圖像分割方法中的“C-V”模型(Active contour without edge)[14]。如圖5和圖7分別是原始的腦部和體部CT圖，圖6和圖8分別是對(duì)應(yīng)分割后的腦部和體部CT圖。從實(shí)驗(yàn)結(jié)果來看，腦部CT圖基本可以將腦室正確地區(qū)分出來，體部CT圖可以非常正確地分割出器官的位置，“C-V”模型用于CT圖像的分割能夠得到較為理想的結(jié)果。

圖5 腦部原始CT圖

圖6 分割后的腦部CT圖

圖7 體部原始CT圖

圖8 分割后的體部CT圖

2 小結(jié)

基于偏微分方程的圖像分割方法確實(shí)是一類很有發(fā)展前景的方法：

（1）使用偏微分方程可以用廣義上連續(xù)的二維函數(shù)對(duì)圖像進(jìn)行建模，從而可對(duì)圖像進(jìn)行求導(dǎo)求積分運(yùn)算，使圖像處理問題規(guī)范化。

（2）可以很好地利用現(xiàn)在一些非常完備的數(shù)值分析和偏微分方程計(jì)算方法來進(jìn)行運(yùn)算[15-16]。

（3）使用偏微分方程可以使圖像分割的速度，準(zhǔn)確性和穩(wěn)定性大幅度提高。

醫(yī)學(xué)圖像分割技術(shù)仍然是制約醫(yī)學(xué)圖像處理中其他相關(guān)技術(shù)發(fā)展和應(yīng)用的一個(gè)瓶頸。醫(yī)學(xué)圖像分割技術(shù)是開展醫(yī)學(xué)圖像在臨床和生命科學(xué)研究領(lǐng)域廣泛應(yīng)用的一個(gè)先行程序，CT圖像分割由于被用于臨床醫(yī)療，因此CT圖像分割的準(zhǔn)確性更為重要?，F(xiàn)在，醫(yī)學(xué)圖像分割的算法大多是把多種理論結(jié)合起來應(yīng)用，以達(dá)到相互補(bǔ)充的目的。目前，新的分割方法的研究主要以自動(dòng)、精確、快速、自適應(yīng)性和魯棒性等幾個(gè)方向作為研究目標(biāo)，偏微分方程方法在這些方面的應(yīng)用具有廣闊的前景，新的偏微分方程模型的研究有待進(jìn)一步努力。

[1] Wells W M,Grimson W E L,Kikins R,et al.Adaptive segmentation of MRI datap[J].IEEE Trans.on Medical Imaging, 1996,15(4):429-442.

[2] 田捷,包尚聯(lián),周明全.醫(yī)學(xué)影像處理與分析[M].北京:電子工業(yè)出版社,2003:35-37.

[3] Dzung L Pham,Chenyang Xu,Jerry L Prince.A Survey of Current Methods in Medical Image Segmentation[R].Technical Report JHU/ECE 99-01,Johns Hopkins Univ,1998.

[4] James S Duncan,Nicholas Ayache.Medical Image Analysis: Progress over Two Decades and the Challenges Ahead[J].IEEE Transaction on patter analysis and machine intelligence,2000,22 (1):181-204.

[5] 李景,等.醫(yī)學(xué)圖像分割技術(shù)[J].生物醫(yī)學(xué)工程學(xué),2006,23(4): 891-894.

[6] 章毓晉.圖像分割[M].北京:科學(xué)出版社,2001.

[7] 章毓晉.圖像處理和分析[M].北京:清華大學(xué)出版社,1999.

[8] 吳月娥,周康源,李傳富,等.基于level sets的醫(yī)學(xué)圖像分割[J].北京生物醫(yī)學(xué)工程,2006,(3):240-243.

[9] Castlemen K R.Digital Image Processing[M].Prentice Hall,1996.

[10] Caselles V,Kimmel R,Sapiro G.Geodesic active contours[J]. International Journal of ComputerVision,1997,22(1):61-79.

[11] Kass M,Witkin A,Terzopoulos D.Snakes:active contour models[J].International Journal of Computer Vision,1988,(33): 321-331.

[12] Chan T F ,Vese L A.Active contours without edges[J].IEEE Trans.Image Processing,2001,10(2):266-277．

[13] 李榮華.偏徽分方程數(shù)值解法[M].北京:高等教育出版杜, 2005.

[14] 王大凱,侯榆青,彭進(jìn)業(yè).圖像處理的偏微分方程方法[M].北京:科學(xué)出版社,2009.

[15] 李玲,等.基于改進(jìn)的Snake模型分割肝臟的CT圖像[J].中國醫(yī)療設(shè)備,2010,25(2):19-20.

[16] 張宜,陳剛.基于偏徽分方程的圖像處理[M].北京:高等教育出版,1990:629-639.

CT Image Segmentation Based on Partial Differential Equation and Programming in MATLAB

ZHOU Ying-jie,WU Hai-jian
The Affiliated Tumor Hospital of Guangzhou Medical University, Guangzhou Guangdong 510095, China

CT images are widely used in medical science. CT image segmentation is important for disease area extraction, tissue location, tissue measuring, and three-dimensional reconstruction. Recently, image segmentations based on partial differential equation(PDE),which is one of the novel and efficient segmentation ways, are gradually turned into research hotspot. In this article, we study in GAC model and C-V model based-PDE image segmentation. The arithmetic of the paper is carried out by programming in matlab,and is verified effectively.

CT images; CT image segmentation; partial differential equations; GAC model; C-V model

R318

B

10.3969/j.issn.1674-1633.2011.05.060

1674-1633(2011)05-0159-03

2010-08-23

2011-03-23

作者郵箱：zhouyingjie@tom.com