利用改進CV模型連續(xù)水平集算法的核磁共振乳腺圖像分割

王芳梅，范虹，Yi WANG

（1.陜西師范大學計算機科學學院，710062，西安；2.康奈爾大學生物醫(yī)學工程系，14853，美國紐約伊薩卡）

隨著X射線鉬靶攝影、超聲成像、X射線計算機斷層攝影（computerized tomography，CT）、核磁共振成像（magnetic resonance imaging，MRI）等影像技術的不斷發(fā)展，乳腺疾病的計算機臨床輔助診斷技術得到了廣泛應用。動態(tài)增強核磁共振成像（Dynamic Contrast-Enhanced Magnetic Resonance Imaging，DCE-MRI）是通過靜脈注射對比劑無損地評價組織和腫瘤血管特性的一種功能性成像方法［1］。MRI因具有無創(chuàng)、多層、多參數、多序列和較高軟組織分辨能力的特性，因而成為診斷乳腺疾病的首選方式之一，MR乳腺圖像分割在乳腺疾病的診斷和量化分析中起著越來越重要的作用。但是，由于MR乳腺圖像存在信息量大、噪聲復雜、邊界模糊、灰度不均勻等現象，使得乳腺活檢圖像的分割工作仍然存在一定的困難。現有的醫(yī)學圖像分割方法主要有閾值法、區(qū)域增長法、特征聚類法和主動輪廓法。

屬于主動輪廓模型范疇的水平集方法［2］因能夠直接獲取閉合分割線、并能有效處理圖像分割過程中的拓撲問題而贏得了關注，目前已得到了大量應用［3-5］。Ball等人利用水平集方法實現了乳房X射線圖像分割，通過利用基于Chan-Vese（CV）模型的離散水平集方法引入目標的區(qū)域特征信息，解決了無明顯邊界信息的目標分割難題，同時解決了傳統CAD系統僅能處理邊緣分析問題的局限性［6］。Lin等人提出了結合高斯金字塔的多分辨率水平集框架，解決了模糊超聲圖像難分割的問題［7］。郭禮華等人利用一種基于高斯矢量場（GVF）的多分辨率水平集圖像分割方法，實現了實際場景的人臉提取，通過引入GVF模型，擴大圖像力的作用范圍，有效解決了弱邊界問題［8］。余瑞星等人利用多分辨率水平計算法，實現了海岸線SAR圖像的分割［9］。

為了進一步研究MR乳腺圖像分割，根據傳統CV模型［10］，本文提出一種基于改進CV模型的連續(xù)水平集算法，利用B樣條基函數［11］將傳統離散水平集函數表示成連續(xù)形式，用解決B樣條空間的變分問題代替水平集函數更新的計算問題，通過在CV模型中引入偏移量構造新的圖像分割能量模型，可以在圖像分割過程中避免水平集函數的重新初始化和分割陷入局部極小值現象，提高分割效率，并有效抑制噪聲，從而獲得準確、穩(wěn)定的MR乳腺圖像分割結果。

1 傳統水平集方法

1.1 離散水平集函數的表示

假設二維平面曲線為C（r，t），其中r為弧長參數，t為時間參數，它對應的三維水平集函數為φ（x，y，t），則曲線可以通過零水平集函數表示為：通過符號距離函數（signed distance function，SDF）表示的水平集函數為

式中：d（x，y）是空間區(qū)域中一點到曲線C的距離，點（x，y）在曲線C的內部時取正號，在外部時取負號。符號距離函數的取值示意圖見圖1。

圖1 符號距離函數取值示意圖

1.2 CV能量函數

在 Mumford-Shah模型［12］的基礎上，Chan和Vese提出了Chan-Vese（CV）模型。該模型是基于區(qū)域信息的圖像分割能量模型，既可以檢測出有梯度信息的目標邊界，也可以檢測出無梯度信息的目標邊界；同時，不僅可以分割出連續(xù)邊界的圖像，還可以分割出不連續(xù)邊界的圖像。用分段常量表示圖像，假設c1、c2為常量，分別表示圖像目標和背景區(qū)域，則CV能量函數的表達式為

式中：μ≥0、ν≥0、λ1＞0、λ2＞0均為常數系數；第1項和第2項分別是長度、面積規(guī)則項；后2項驅動閉合曲線項的目標邊界演化。u0：Ω→R是原圖像區(qū)域，初始輪廓線用C表示，c1、c2是通過給定閉合曲線C內、外區(qū)域的灰度均值表示的。

2 基于改進CV模型的連續(xù)水平集算法的實現

2.1 連續(xù)水平集函數的表示

B樣條是一種簡單的基函數，整數結點處的函數值可以表示該點處基函數的系數。水平集函數可用B樣條基函數的線性組合形式表示為

式中：βn（·）是n階B樣條基函數；h為空間尺度；B樣條的結點是由圖像定義的空間域Ω中的網格結點；c（i，j）表示 B樣條基函數的系數集。

由于水平集函數是由B樣條基函數的線性組合表示的，因此，水平集函數不斷更新的過程可以用簡單的卷積操作實現。水平集的每一步演化都可以視為經過B樣條核函數濾波作用的結果，該濾波器具有較高的抗噪性，可保證在圖像分割中對噪聲的抑制作用。

2.2 基于改進CV模型的圖像分割模型

Chan和Vese在Mumford-Shah能量式基礎上改進的圖像分割方法——CV模型結合了圖像的全局灰度統計信息，通過極小化水平集函數的能量函數來達到圖像分割的目的。為了避免在求解函數極小值過程中出現局部極小值的現象，Chan和Vese引入了Heaviside函數，但是卻存在難以設置曲線演化停止條件的困難。

為了實現乳腺MR圖像的分割，本文對CV模型中的擬合項引入偏移量得到α-CV模型，并將其作為能量函數。將此能量函數用于分割實驗，可以去除局部極小值現象并設置有效的曲線演化停止條件。能量函數的表達式為

為了求式（5）的極小值，需利用關于水平集函數的梯度下降流方程式

式中：α是任意小的非負常量；Ω是圖像u0的定義域。關于水平集函數φ的能量函數式的極小化得到的零水平集將目標區(qū)域與背景區(qū)域分割成2部分。通過式（5）等號右側第2項的負號可見，該能量泛函明顯區(qū)別于CV模型。通過擴大水平集函數φ的比例，使得能量泛函值可以不依賴于φ值符號的改變而變化，從而避免了能量泛函值陷入局部極小值現象。式中修改Heaviside函數為位移Heaviside函數，用于控制決定能量泛函值變化的水平集函數值的范圍。當唯一變量α=0時，能量項中的第1項僅在水平集函數為正值的情況下才發(fā)生變化，而能量項中的第2項僅在水平集函數為負值的情況下才發(fā)生變化。即使u0（X）≠c1，u0（X）≠c2，當求能量項的極小值時，會使得水平集函數的值趨向于0，從而發(fā)生曲線局部收斂。為了避免這種現象，本文引入帶有位移量α的Heaviside函數，具體作用如下：

（1）控制水平集函數的演化，根據當前水平集函數確定能量函數的極小化，從而達到全局極小值的目的；

（2）制約水平集函數的取值，從而保證求函數極小值的過程能得到穩(wěn)定解。

當u0（X）≠c1、u0（X）≠c2時，水平集函數應該趨于某個非零值，這樣零水平集才能實現圖像域中目標與背景的分離。為了達到這一目的，將Heaviside函數按照變量α進行位移。對式（6）僅保留第1項，對能量泛函極小化，當任意時刻φ＞-α時，φ值趨于-α；對式（6）僅保留第2項，對能量泛函極小化，當任意時刻φ＜α時，φ值趨于α；當這2項全部保留時，在水平集函數值滿足-α≤φ≤α的圖像域中，由｜u0-c1｜2、｜u0-c2｜2這2項的取值決定水平集函數φ是否趨于α或-α，由此達到期望的分割目的。顯然，通過式（5）可得：當能量泛函被極小化時，任意處于｜u0-c1｜2＜｜u0-c2｜2圖像域中的像素點都被劃分到區(qū)域r：φ（r）=-α中；相反，任意處于｜u0-c1｜2＞｜u0-c2｜2圖像域中的像素點都被劃分到區(qū)域r：φ（r）=α中。

假設圖像定義域為2相，且滿足R1∈Ω：｜φ（R1）｜＞α和R2∈Ω：｜φ（R2）｜≤α。為了方便理解，取λ1=λ2=1，定義能量泛函積分的表達式為

可見式中φH（α+φ）和φH（α-φ）是由水平集控制的，即｜u0-c1｜2、｜u0-c2｜2僅依賴于零水平集。當在R1中時，有Γ（φ）=｜u0-c1｜2φ，Γ（φ）隨著φ的減小而減小，直到水平集到達目標邊界，即φ=α時，Γ（φ）達到穩(wěn)定狀態(tài)；相反，當在R2中時，Γ（φ）隨著φ的減小而減小，直到水平集到達目標邊界，即φ=-α時，Γ（φ）達到穩(wěn)定狀態(tài)。因此，在式（6）的作用下，所有滿足｜φ｜＞α區(qū)域中的水平集隨著｜φ｜的減小而演化一次，并且控制水平集函數處于-α≤φ≤α區(qū)域。

3 實驗結果與分析

3.1 仿真數據實驗結果與分析

為了驗證算法的有效性，在Matlab（R2011b）編程環(huán)境下進行仿真實驗，實驗平臺為 Windows XP，CPU為Pentium 4，主頻為3.0GHz，內存為1GB。圖像大小為256×256像素，灰度級為256。迭代步長為3。為了進一步簡化正則化函數δ（·），參數ε取1。

分割結果如圖2所示：圖中第1列是原始圖像，第2列是加入均值和方差分別為0和0.01的高斯噪聲后的圖像（稱為含噪圖像Ⅰ），第3列是加入均值和方差分別為0和0.1的高斯噪聲后的圖像（稱為含噪圖像Ⅱ）；第1行是分割前的圖像，第2行是使用傳統基于CV模型的離散水平集算法（以下簡稱傳統算法）分割后的圖像（ν=255×255，迭代100次），第3行是在h=2、ν=0、迭代3次條件下使用本文提出的基于改進CV模型的連續(xù)水平集算法（以下簡稱本文算法）的分割結果，第4行是在h=3、ν=0、迭代3次條件下采用本文算法的分割結果。

圖2 數字圖像分割對比圖

從抗噪性能方面看，由于傳統離散水平集算法對曲率的高度依賴，使得平滑噪聲的效果存在一定的限制，信噪比越低的圖像平滑效果越差，對噪聲敏感性較強。本文算法可以有效降低對噪聲的敏感性，和傳統離散水平集方法相比，即使分割較低信噪比的圖像，也能得到較好的分割效果。由圖2第3行和第4行的分割結果可以看出，文中所選基函數參數h的取值對分割效果的影響呈正相關，即在相同條件下，h取值越大對噪聲的平滑效果越好，得到的圖像分割效果越理想。

從分割精度方面看，由圖2第2、第3和第4行相應列的分割結果可以看出：傳統CV模型存在容易陷入局部極小值的現象，得不到全局圖像分割結果；本文采用的改進CV模型擬合項的能量函數可以避免出現局部極小值現象，得到較好的全局圖像分割結果，并且在同樣的條件下達到有效收斂，完成曲線演化，得到全局目標輪廓線。

3.2 臨床數據實驗結果與分析

為了實現對MR乳腺圖像的分割，本節(jié)將對4個病人的臨床實驗數據（13個）用基于CV模型的離散水平集算法和本文算法進行分割對比，通過仿真實驗證明本文算法分割圖像的有效性。

實驗中所用的數據是由康奈爾大學醫(yī)學成像實驗室提供、并且由病人授權本論文使用的乳腺臨床MRI數據，圖像是通過使用T1加權獲得的帶有飽和脂肪梯度回波力的矢狀面切片序列，采樣間隔時間為8.1ms，回波時間為4ms，翻轉角度為30°。成像所用儀器為德國埃朗根西門子公司的具有專用雙乳房線圈的1.5T全身核磁共振成像系統，病人在成像時處于俯臥姿勢。

圖3～圖5是針對3位病人在相同采樣間隔時間下獲取的冠狀面、矢狀面成像的12例乳腺MR圖像，采用本文基于改進CV模型的連續(xù)水平集算法和傳統的基于CV模型的離散水平集算法的分割結果對比，圖像大小為512×512像素，平均灰度級分別為1 125、2 215和667。水平集函數參數h=3，ν=0，迭代10次。

圖3 矢狀面乳腺MR圖像A的分割對比圖

圖4 矢狀面乳腺MR圖像B的分割對比圖

圖5 冠狀面乳腺MR圖像的分割對比圖

圖6是用本文算法和傳統算法對橫斷面成像圖像的分割結果對比圖，圖像大小為256×256像素，灰度級為499。圖中第1列是原始圖像，第2列是加入均值和方差分別為0和0.01的高斯噪聲后的圖像（即含噪圖像Ⅰ），第3列是加入均值和方差分別為0和0.1的高斯噪聲后的圖像（即含噪圖像Ⅱ）；第1行是分割前的圖像，第2行是傳統算法的分割結果（ν=499×499，迭代200次）；第3行和第4行是本文算法的分割結果，其中第3行參數為h=2，ν=0，迭代5次；第4行調整空間尺度h為3，即參數為h=3，ν=0，迭代5次。

圖6 橫斷面乳腺MR圖像的分割對比圖

根據實驗結果可以看出，通過B樣條基函數的連續(xù)水平集函數表示，大大降低了迭代次數，可通過靈活調整基函數的空間尺度參數h達到抑噪效果，且傳統算法的局部極小值現象在本文算法中已經解決，即通過本文算法可以得到全局圖像分割結果。由圖3～圖5可以看出，本文算法比傳統算法獲取的分割結果更為光滑且接近目標輪廓。比較圖6第3行和第4行可以證明，本文利用B樣條基函數表示水平集函數時，針對特定圖像，可以通過適當調整空間尺度h的值得到最佳分割效果；比較第2行和第3行的最后一列分割圖可以看出，本文算法對較高信噪比圖像的分割效果明顯好于傳統算法的分割效果，即具有更好的抗噪性。

4 結 論

本文針對MR乳腺圖像的分割難題進行了研究，提出一種基于改進CV模型的連續(xù)水平集算法，采用B樣條基函數表示連續(xù)水平集函數，相比傳統離散水平集函數大大降低了計算過程的復雜性，避免了水平集函數的重新初始化；基于改進CV模型的圖像分割能量模型避免了分割陷入局部極小值現象，并且可有效抑制噪聲。通過實驗對仿真數據和臨床數據進行了分割，結果表明：本文提出的基于改進CV模型的連續(xù)水平集算法可以準確、穩(wěn)定地實現低信噪比、弱邊界MR乳腺圖像的分割，并具有一定的抗噪性能。

［1］ KNOPP M V，GIESEl F L，MARCOS H，et al.Dynamic contrast-enhanced magnetic resonance imaging in oncology［J］.Magn Reson Imaging，2001，12（4）：301-308.

［2］ OSHER S，SETHIAN J A.Fronts propagating with curvature dependent speed：algorithms based on Hamilton-Jacobi formulations［J］.Journal of Computational Physics，1988，79（1）：12-49.

［3］ 龔永義，羅笑南，黃輝，等.基于單水平集的多目標輪廓提取 ［J］.計算機學報，2007，30（1）：121-128.GONG Yongyi，LUO Xiaonan，HUANG Hui，et al.Multi-objects extracted based on single level set［J］.Chinese Journal of Computers，2007，30（1）：121-128.

［4］ 王曉飛，龐全.基于圓形約束快速水平集的原生質體細胞分割 ［J］.中國圖象圖形學報，2013，18（1）：55-61.WANG Xiaofei，PANG Quan.Protoplasm somatic cells segmentation based on circle dependent fast levelset segmentation［J］.Journal of Image and Graphics，2013，18（1）：55-61.

［5］ 潘改，高立群，趙爽.基于局部熵的主動輪廓模型［J］.中國圖象圖形學報，2013，18（1）：78-85.PAN Gai，GAO Liqun，ZHAO Shuang.Active contour model driven by local entropy energy［J］.Journal of Image and Graphics，2013，18（1）：78-85.

［6］ BALL J E，BRUCE L M.Level set-based core segmentation of mammographic masses facilitating three stage（core，periphery，spiculation）analysis［C］∥Proceedings of the 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society，2007.Piscataway，NJ，USA：IEEE，2007：819-824.

［7］ LIN Ning，YU Weichuan，DUNCAN J S.Combinative multi-scale level set framework for echocardiographic image segmentation［J］.Medical Image Analysis，2003，7（4）：529-537.

［8］ 郭禮華，李建華，袁小彤，等.基于高斯矢量場的多尺度水平集圖像分割算法 ［J］.上海交通大學學報，2005，39（8）：129-133.GUO Li-hua，LI Jian-hua，YUAN Xiao-tong，et al.The multi-scale image segmentation algorithm based on Gaussian vector field（GVF）snake model［J］.Journal of Shanghai Jicaotong University，2005，39（8）：129-133.

［9］ 余瑞星，李嚴俊，張科.基于小波變換的多尺度水平集算法研究 ［J］.光子學報，2007，36（2）：372-375.YU Rui-xing，LI Yan-jun，ZHANG Ke.Multiresolution level set studying based on wavelet transform ［J］.Acta Photonica Sinica，2007，36（2）：372-375.

［10］CHAN T，VESE L.Active contours without edges［J］.IEEE Transactions on Image Processing，2001，10（2）：266-277.

［11］OLIVIER B，DENIS F，PHILIPPE T，et al.Variational B-spline level-set：a linear filtering approach for fast deformable model evolution ［J］.IEEE Transactions on Image Processing，2009，18（6）：1179-1191.

［12］MUMFORD D，SHAH J.Optimal approximations by piecewise smooth functions and associated variational problems［J］.Communications on Pure and Applied Mathematics，1989，42（7）：577-685.