B超圖像在醫學臨床診斷中有著重要應用,但廣泛存在的灰度分布不均勻、對比度低、偽影和噪聲干擾以及目標邊界模糊等問題,給自動分割帶來了困難。本文在區域水平集模型的基礎上定義反映輪廓線像素點對目標/背景兩個區域隸屬度的因子,通過概率分布估計模型計算和比較各像素點的隸屬度,以此為依據對像素點進行區域歸屬判別,由區域水平集迭代獲得連續光滑的曲線。本文將B超圖像目標分割看作對感興趣目標區域的局部分割,將水平集的計算求解約束到局部范圍,從而減少計算量。實驗結果表明與幾種水平集模型相比,本文方法對所測試的B超圖像的分割在精度和速度上均有一定的改進。
引用本文： 楊誼, 喻德曠, 申洪. 隸屬度區域水平集方法提取B超圖像病灶. 生物醫學工程學雜志, 2015, 32(4): 779-787. doi: 10.7507/1001-5515.20150142 復制
引言
超聲診斷是將超聲檢測技術應用于人體，通過測量了解生理或組織結構的數據和形態，從而發現疾病的一種診斷方法，主要用于腫物、畸形、結石及其他能引致局部結構發生明顯形態改變的疾病，具有無創、方便、廉價、實時等優點。其中應用最廣的B型超聲(B-type ultrasound,BU)檢查對液體與實質性組織具有顯著的圖像差別，通過黑白之間的不同灰度來反映不同組織回聲的強弱，對囊性液體與實質性病變組織具有良好的鑒別能力[1]。提取病灶區域是醫學診斷和制定治療計劃的先行步驟。醫學圖像分割根據區域間的相似或不同，把圖像分割成若干區域，提取感興趣組織或部位。由于超聲成像原理、被測生物體組織對聲阻抗分布不勻及在檢測過程中可能發生的運動等因素影響，B超圖像經常存在灰度不均勻、對比度較低、目標邊界模糊、偽影和噪聲干擾等現象，給自動分割帶來了較大的困難。
1 超聲分割方法的探索
1.1 分割方法概述
針對這些問題，近年來研究人員做了大量工作，所使用的方法主要有區域增長法、動態閾值分割法、模糊聚類法、分水嶺法、隱馬爾可夫模型(hidden Markov models，HMM)融合法、最大估值算法(estimation maximization，EM)和高斯混合分布模型、小波變換法等，提出了針對一些具體問題的解決方案。基于變分的水平集方法是近年來研究頗為活躍的一個分支，它將圖像分割問題表達為能量函數的最小化問題，由變分原理將其轉化為偏微分方程的求解[2]。相比于傳統的區域分割方法，變分方法可以在能量函數中加入與圖像內容有關的約束條件等綜合信息，獲得更加自然準確的分割效果。
1.2 水平集方法的應用評述
水平集方法的主要思想是將移動的界面作為零水平集嵌入高一維的水平集函數中，由封閉超曲面零水平集的演化方程得到水平集函數的演化方程，通過水平集函數曲面的演化來隱含地求解輪廓曲線，結果得到目標輪廓閉合曲線。Chan-Vese(CV)模型[3]是區域模型的代表，設圖像I(x,y)的定義域為Ω，設圖像I(x,y)被輪廓曲線C劃分為目標(曲線C的內部)和背景(曲線C的外部)兩個同質區域，平均灰度分別為Cin 和Cout，設計基于整個圖像的區域能量函數，通過面積和長度約束能量項保持曲線的連續性和平滑性，能量泛函為：
| ${E_{CV}}\left( {{c_{{\rm{in}}}},{c_{{\rm{out}}}},C} \right) = {E_{{\rm{global}}}} + {E_{{\rm{length}}}} + {E_{{\rm{area}}}} = {\lambda _1}\int\limits_{{\Omega _{{\rm{in}}}}} {{{\left| {I\left( {x,y} \right) - {c_{{\rm{in}}}}} \right|}^{\rm{2}}}{\rm{d}}x{\rm{d}}y} + {\lambda _2}\int\limits_{{\Omega _{{\rm{in}}}}} {{{\left| {I\left( {x,y} \right) - {c_{{\rm{in}}}}} \right|}^{\rm{2}}}{\rm{d}}x{\rm{d}}y} + \mu \left| C \right| + \gamma A$ |
CV模型利用圖像的全局灰度信息而非邊緣信息，能夠有效分割出模糊或部分缺失的邊緣；當兩種區域的灰度差別較大時具有較強的抗噪能力；對初始輪廓線位置不敏感；經嚴格證明能夠收斂到全局最優。但是，它的前提條件是圖像灰度分布均勻、前景背景有明顯差異，而B超圖像普遍對比度不高，分割結果不夠準確。另外，μ的取值為經驗式的，取值小時能較準確刻畫病灶不規則邊界，但輪廓線常產生大量零碎冗余，取值大時輪廓線形狀光滑卻容易偏離目標的不規則邊界。
近年來針對水平集在醫學中的運用出現了一些較好的改進方案。如文獻[4]采用了圖像平滑和圖像分割同時進行的策略，可以自動分割出多種類型的區域，能自動處理拓撲變化，但由于平滑和分割同時進行導致計算量較大，僅適用于圖像內容比較簡單的情況。文獻[5]設計了一種區域可擴展自適應能量函數(region-scalable fitting energy，RSF)以構造水平集模型，對局部區域內輪廓線的擬合能量函數捕獲圖像的局部區域信息，能夠較好地處理灰度不勻的情況，但運算代價比CV模型顯著增加，使用的核函數的方差為定值，沒有區別考慮不同區域的強度變化，并且容易陷入局部極值解。文獻[6]提出可選擇的局部或全局能量函數(selective local or global segmentation，SLGS)，運用可選擇的二值高斯過濾糾正分割方案，在弱邊界處逼停演化的曲線，具有較好的抗噪性，并且能夠較合理地提取弱邊界，但采用局部區域均值描述區域灰度仍存在一定誤差，這種情況在低對比度B超圖像中經常出現。文獻[7]提出基于局部高斯分布的自適應能量函數(local Gaussian distribution fitting energy，LGF)，將高斯核函數作為權值來刻畫局部區域灰度變化情況，能夠較準確地表現局部的不規則變化，但精度由核函數的半徑確定，不可改變。文獻[8]采用局部統計信息引導CV模型演化，思路與文獻[7]異曲同工，其優點是采用了快速水平集方法，建立存儲曲線內外鄰點的鏈表，根據曲線上點的鄰域的適應能量場對最佳逼近目標邊界位置的收斂條件進行定義，實現演化曲線只需更新這兩個表，計算代價減小了，但維系鏈表比使用數組代碼復雜度增大，精度仍然受到局部區域取值的限制。
本文針對B超圖像特點，利用區域水平集模型的優點，提出一種新的水平集方法，定義反映像素點對目標/背景兩個區域隸屬度的因子，通過概率分布估計模型計算和比較像素點的隸屬度，以此為依據對像素點進行區域劃分，以提高水平集曲線定位B超圖像模糊邊界的能力；同時本文將B超圖像目標分割看作對感興趣目標區域的局部分割，將水平集的計算求解約束到局部分割范圍，減少計算量，提高分割速度，實現病灶等特定目標區域的高效提取。
2 隸屬度區域水平集模型算法
2.1 本文高斯混合模型和水平集模型結合的思路
高斯分布(Gaussian distribution)描述了一種圍繞某個單值聚集分布的隨機變量，是統計學中最廣泛應用的一類分布。中心極限定理表明，采樣的均值近似服從高斯分布，高斯分布的信息熵在所有的已知均值及方差的連續分布中最大，這使得它成為一種均值以及方差已知的分布的自然選擇。高斯分布是一個單極值分布，不能對多極值復雜分布提供一個較好的近似，為此提出了高斯混合模型(Gaussian mixture model，GMM)，把數據看作是從數個高斯分布中生成出來的，將一個事物分解為若干的基于高斯概率密度函數(正態分布曲線)形成的模型，以得到更精確的量化模型。GMM的概率密度函數就是幾個高斯概率密度函數的加權和。圖像灰度直方圖反映的是圖像中某個灰度值出現的頻次，可以認為是圖像灰度概率密度的估計。如果圖像所包含的目標區域和背景區域在灰度上有明顯差異，那么該圖像的灰度直方圖呈現雙峰形狀，分別對應目標和背景的中心灰度。而復雜的圖像如醫學圖像、遙感圖像等，直方圖一般是多峰的。通過將直方圖的多峰特性看作是多個高斯分布的疊加，可以解決圖像的分割問題。
本文將GMM和水平集模型結合起來，定義待處理B超圖像為I(x)=a(x)r(x)，其中r(x)為理想圖像(灰度均勻分布、無噪、無偽影)，I(x)為實際待處理圖像，a(x)為偏差因子。假設理想圖像r(x)中像素點的強度分布滿足GMM：
| $\pi \left( x \right) = \sum\limits_{i = 1}^M {{\omega _i}g\left( {x;{m_i},{\sigma _i}} \right)} $ |
其中M為分割結果中的類別(區域)數目，本文取M=2(目標和背景)，ωi為先驗概率，滿足，g(x;mi,σi)為均值為mi、協方差為σi的Gaussian分布：g(x;mi,Σi)=。
式(3)表示的概率模型似然函數為
| $L\left( {u\left( x \right)} \right) = \prod\limits_{i = 1}^M {\prod\limits_{{\Omega _i}} {{\omega _i}{p_i}\left( {u\left( x \right)} \right)} } $ |
最大似然法(maximum likelihood,ML)是一種點估計法，當從模型總體隨機抽取n組樣本觀測值后，最合理的參數估計量應該使得該n組樣本觀測值的概率最大。根據最大似然估計法，求總體參數u(x)的極大似然估計值的問題就是求似然函數L(u(x))的最大值(或-L(u(x))的最小值)問題，可通過解=0來解決。因為logL是L的增函數，兩者在u(x)的同一值處取得最大值，就把問題轉化為logL的極值求解。
2.2 本文隸屬度區域水平集方法
設C={Cin,Cout}是圖像I關于輪廓線C的兩個劃分，Cin代表輪廓線內部即目標區域，Cout代表輪廓線外部即背景區域，根據文獻[9]，設p(Cin)和p(Cout)分別為目標區域和背景區域關于像素點x的先驗概率，則有p(Cin)=p(Cout)。像素點x屬于目標或背景的概率分別為
| $\begin{array}{l} p\left( {C\left. {_{{\rm{in}}}} \right|x} \right) = \frac{{p\left( {{C_{{\rm{in}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right)}}{{p\left( {{C_{{\rm{in}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right) + p\left( {{C_{{\rm{out}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right)}}\\ p\left( {C\left. {_{{\rm{out}}}} \right|x} \right) = \frac{{p\left( {{C_{{\rm{out}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right)}}{{p\left( {{C_{{\rm{in}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right) + p\left( {{C_{{\rm{out}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right)}} \end{array}$ |
其中，u(x)為I(x)中像素點x的灰度值，定義關于每個像素點x的區域隸屬度因子p(u(x)｜Cin)和p(u(x)｜Cout)，分別表示在目標類與背景類條件下出現u(x)的概率，它們服從互相獨立的高斯分布。由p(Cin)p(u(x)｜Cin)+p(Cout)p(u(x)｜Cout)為常數，p(Cin)=p(Cout)=0.5，將式(5)簡化為：
| $\begin{array}{l} p\left( {C\left. {_{{\rm{in}}}} \right|x} \right) = {\lambda _1}p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right) = \frac{\lambda }{{{\sigma _1}}}\exp \left\{ { - \frac{{{{\left( {x - {m_1}} \right)}^2}}}{{2\sigma _1^2}}} \right\}\\ p\left( {C\left. {_{{\rm{out}}}} \right|x} \right) = {\lambda _2}p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right) = \frac{\lambda }{{{\sigma _2}}}\exp \left\{ { - \frac{{{{\left( {x - {m_2}} \right)}^2}}}{{2\sigma _2^2}}} \right\} \end{array}$ |
其中。
根據區域隸屬度因子進行像素點的區域歸屬判定，當p(Cin｜x)>p(Cout｜x)時，像素點x劃歸為目標類；p(Cin｜x)<p(Cout｜x)時，像素點x劃歸為背景類；p(Cin｜x)=p(Cout｜x)時，則像素點x為目標類和背景類的分界點。將目標和背景區域的符號取相反(內負外正)以實現求差運算，定義能量函數為：
| $\begin{array}{l} E\left( C \right) = - \lambda \int\limits_\Omega {p\left( {C\left. {_{{\rm{in}}}} \right|x} \right){\rm{d}}{R_{{\rm{in}}}}} + \lambda \int\limits_\Omega {p\left( {C\left. {_{{\rm{out}}}} \right|x} \right){\rm{d}}{R_{{\rm{out}}}}} + \mu \left| C \right| + \gamma A = - \lambda \int\limits_\Omega {\left[ {\log \frac{{{\sigma _1}}}{{{\sigma _2}}} - \frac{{{{\left( {u\left( x \right) - {m_1}} \right)}^2}}}{{\sigma _1^2}} + \frac{{{{\left( {u\left( x \right) - {m_2}} \right)}^2}}}{{\sigma _2^2}}} \right]} {\rm{d}}x + \\ \mu \left| C \right| + \gamma A \end{array}$ |
顯然，當p(Cin｜x)=p(Cout｜x)時式(6)達到最小值，水平集停止在目標和背景的邊界上，得到分割輪廓線。使用高斯分布計算概率分布值需要先固定分割區域，求出各區域內所有像素灰度的均值mi和方差σi作為概率分布參數。在水平集曲線演化過程中分割區域是隨輪廓線時刻變化的，每迭代一次就要根據當前的分割情況采樣所有點，計算不同區域的概率分布參數計算量很大。考慮到B超圖像中不同的組織具有不同的灰度值，每個區域內部也存在強度不均勻性，以小鄰域信息代替全局信息能夠更準確地刻畫圖像局部特征。在每個像素點x的小鄰域內計算圖像的統計信息，以更準確地反映出這種變化的強度分布。引入以像素點x為中心，大小為k×k的窗口函數wk(x)，把窗口wk(x)中曲線C的內部區域和外部區域的灰度均值分別記為minw(x)和moutw(x)，方差分別記為σinw2(x)和σoutw2(x)，內外兩個區域的像素點個數分別為ninw(x)和noutw(x)，灰度均值和方差的計算式為：
| $\begin{array}{l} {m_{{\rm{inw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{in}}}}\left( x \right)} {I\left( x \right)} }}{{{n_{{\rm{inw}}}}}},{m_{{\rm{outw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{out}}}}\left( x \right)} {I\left( x \right)} }}{{{n_{{\rm{outw}}}}}},\\ {\sigma _{{\rm{inw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{in}}}}\left( x \right)} {{{\left( {I\left( x \right) - {m_{{\rm{inw}}}}\left( x \right)} \right)}^2}} }}{{{n_{{\rm{inw}}}}}},\\ {\sigma _{{\rm{outw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{out}}}}\left( x \right)} {{{\left( {I\left( x \right) - {m_{{\rm{outw}}}}\left( x \right)} \right)}^2}} }}{{{n_{{\rm{outw}}}}}} \end{array}$ |
引入水平集表示法，得到隸屬度區域水平集能量演化方程：
| $\frac{{\partial \Phi }}{{\partial t}}\delta \left( \Phi \right)\left[ {\log \frac{{{\sigma _{{\rm{inw}}}}}}{{{\sigma _{{\rm{outw}}}}}} - \frac{{{{\left( {u\left( x \right) - H\left( {\Phi \left( {u\left( x \right)} \right)} \right) * {m_{inw}}} \right)}^2}}}{{\sigma _{{\rm{inw}}}^2}} + \frac{{{{\left( {u\left( x \right) - \left( {1 - H\left( {\Phi \left( {u\left( x \right)} \right)} \right)} \right) * {m_{{\rm{outw}}}}} \right)}^2}}}{{\sigma _{{\rm{outw}}}^2}} + \mu {\rm{div}}\left( {\frac{{\nabla \Phi }}{{\left| {\nabla \Phi } \right|}}} \right) + \gamma } \right]$ |
其中Hε(x)為Heaviside函數，記做Hε(x)=，δε(ф)=H′ε(x)=。這樣就把圖像目標分割問題轉化為極小化能量函數(8)的計算。
3 計算方法的具體實現
3.1 初始輪廓的提取
大多數區域水平集模型是全局分割模型，將圖像中的所有目標分割出來，計算量很大，對于目標針對性強的B超圖像無此必要。本文水平集模型對初始輪廓的設置比較敏感，初始輪廓的設置得當，有助于提高分割準確度和效率。超聲病灶在不同圖像中位置并不固定，如果將初始輪廓統一設置在固定位置、統一形狀(如橢圓或矩形)，會導致初始輪廓與目標存在較大偏離，造成計算時間過長，甚至產生錯誤的分割結果。
經過比較測試，本文采用計算簡單的最大類間方差法(Otsu算法)自動確定圖像分割的初始條件。最佳閾值法按圖像的灰度特性將圖像分成背景和前景兩部分，背景和前景之間的類間方差越大，說明構成圖像的兩部分的差別越大，當部分前景錯分為背景或部分背景錯分為前景時都會導致兩部分差別變小，設灰度圖像灰度級是L，灰度范圍為[0,L-1]，最佳閾值計算公式為：
| $t = \max \left( {{w_b}\left( t \right) * {{\left( {{u_b}\left( t \right) - \bar u} \right)}^2} + {w_f}\left( t \right) * {{\left( {{u_f}\left( t \right) - \bar u} \right)}^2}} \right)$ |
其中，wb為背景比例，ub為背景均值，wf為前景比例，uf為前景均值，
Otsu二值化結果中常存在著許多碎片，需要進行后續處理，本文將數學形態學的開閉運算結合起來對二值圖像進行噪聲濾除，以開運算為主，先腐蝕后膨脹，以消除孤立噪音等小物體、斷裂背景與目標的纖細連接處、平滑較大物體的邊界；以閉運算為輔，先膨脹后腐蝕，填充目標內部細小空洞、連接鄰近物體、平滑邊界。在此基礎上運用八鄰域跟蹤算法勾勒病灶邊界，提取初始輪廓，作為水平集的起點進行演化運算。
圖 1為提取初始輪廓的過程。需要注意的是，由于部分B超圖像內容復雜，當出現多個閉合曲線時需要人工選擇目標區域；當閉合曲線與背景有錯誤粘連時，需要人工截斷錯處的連接點，設置正確的起點，重新進行初始輪廓點的搜索。
圖1
確定水平集初始輪廓的過程
(a)待處理素材；(b)Otsu二值化結果；(c)數學形態學開閉運算的濾波結果；(d)八鄰域跟蹤算法提取的閉合曲線
Figure1. Procedure of getting initial contour for level set(a) image to be processed； (b) result of Otsu binary process； (c) result of mathematical morphology filting； (d) closed curve achieved by eight-neighborhood tracking
3.2 快速高效的局部計算方案
為提高計算效率，本文水平集模型采用局部計算方案，定義基于局部近似符號距離函數，放寬對初始水平集函數的要求，不再要求初始水平集必須是符號距離函數。采用曲線上每個像素點周圍的k×k個鄰居點的聯合概率分布代替該點的概率分布，以減少由于灰度分布不均、對比度低和噪聲干擾帶來的計算誤差，比采用全局信息更能準確地定位模糊邊界。但如果直接對輪廓線上的每個像素點都做對數及求和運算，求解水平集式需要離散化迭代實現求解偏微分方程，計算量較大，曲線演化速度較慢。因此采用局部計算思路以減少運算量和運算時間。文獻[9-10]提出了改進的計算方案以及免除反復初始化水平集的做法，能夠有效加快演化速度，但某些實現技術較為復雜。本文借鑒文獻[11]，采用改進的FTC(Fast two-cycle)算法，在不用求解偏微分方程的基礎上近似實現水平集函數的進化，更新局部近似符號距離函數的定義為
| ${\Phi _0}\left( {x,y} \right) = \left\{ {\begin{array}{*{20}{l}} {k,點x為外部點}\\ {1,點x為輪廓線外領域點}\\ {0,點x在輪廓線上}\\ { - 1,x為輪廓線內領域點}\\ { - k,點x為內部點} \end{array}} \right.$ |
水平集迭代過程是：
(1)如果輪廓線向外運動一個像素，則判斷Lout中的每個點x的能量變化ΔELout(新能量值－原能量值)的符號情況：
①如果ΔELout≥0則不處理；
②如果ΔELout<0，則輪廓線在該點向外擴展一個像素，并且掃描Lin，如果存在x的鄰接點y，則修改ф(y)為-k，即當曲線向外擴展時將原來的內鄰域點修改為內部點。
(2)如果輪廓線向內運動一個像素，則判斷Lin中的每個點x的能量變化 ΔELin(新能量值-原能量值)的符號情況：
①如果ΔELin≥0則不處理；
②如果ΔELin<0，則輪廓線在該點向內擴展一個像素，并且掃描Lout，如果存在x的鄰接點y，則修改ф(y)為k，即當曲線向內擴展時將原來的外鄰域點修改為外部點。
一般情況下，k值取得較大，則分割較為精確，但考慮到計算代價，k值一般取為3或5。每次迭代結束后所得目標輪廓作為下一次迭代的初始輪廓，對圖像進行連續分割，當能量函數的變化小于一定的閾值時，迭代終止，得到分割結果。只需對輪廓線上的點而不是整幅圖像的點進行計算，通過計算能量函數的增減變化來確定輪廓線的擴展或收縮。本文方法的局限是需要保持初始輪廓接近目標邊界，否則就必須擴大鄰域窗口半徑來保證輪廓線形變的正確性，而加大鄰域窗口則意味著運算量的增加。
4 結果
分割實驗待測圖像來自南方醫科大學附屬醫院影像庫，通過對相同B超圖像樣本進行分割來檢驗本文模型的效果，并與幾種性能優良的水平集模型進行比較。編程平臺為Matlab7，計算機配置為Intel? CoreTM(2) Duo CPU，3.0 GHz，2 GB RAM，Windows XP 操作系統。所有實驗采用相同的終止閾值ΔE=10-4或最大迭代次數2 000(考慮到用戶等待時間和實用性)。圖 2~圖 4分別為臨床B超待分割圖像、本文模型初始輪廓、CV模型[3]、RSF模型[5]、SLGS模型[6]、LGF模型[7]、本文模型分割結果和專業醫師手工分割結果。所有模型均按原文獻方法進行相應參數設置，其中公共參數：時間步長Δt=0.1，空間步長Δh=1，對Heaviside函數和Dirac函數進行近似時正則化參數ε=1，長度懲罰項參數μ、鄰域半徑等則根據圖像內容調整為最合適值，以盡可能達到每種模型最優的分割效果。
從分割結果的視覺效果和統計信息兩個方面來比較5種模型在處理B超圖像方面的適用性。
4.1 實驗結果的定性分析
首先看分割結果的視覺效果。圖 2第1行左一為待測樣本411*290像素腎盂積水B超截圖，分割目標為中部位置的暗色積液區域。整幅圖像存在明顯的灰度分布不勻，有多處偽影存在，左側存在大面積偽影，上部和右側存在帶狀偽影，存在散布的噪聲雜質，積液區域的上部有兩處凸起，邊界較為模糊，左側和右側與偽影的交界比較模糊。圖 2第1行左二中的紅色閉合曲線為初始輪廓，并未將病灶完全包括在內，而是與積液區域交疊。以手工分割結果為標準對各個模型的分割結果進行評判：CV模型分割結果較好地處理了灰度不勻的現象，并具有一定的抗噪聲能力；RSF模型分割結果顯示，曲線較容易越過邊界向病灶內部運動，最終收斂到噪聲處，無法準確貼近病灶邊界；SLGS模型分割結果對不規則處的處理較好，但容易越過病灶邊界深入到病灶內部，使得病灶有效提取區比實際略小，同時也把圖像中無關的部分中與病灶灰度分布相似的成分也分割出來，造成不必要的干擾；LGF模型分割結果雖能基本準確地提取病灶邊界，但產生了大量無用的分割碎片，分割結果必須經過后期處理才能用于輔助診斷；而本文方法分割結果在病灶提取方面接近CV模型，對不規則凸起和凹陷處的提取則更為合理，在刻畫模糊邊界處也更為精細準確，能夠反映病灶形態的變化，最接近手工分割結果，但提取的病灶區比實際略大。
圖2
腎盂積水B超圖像分割結果
Figure2.
Segmentation of BU image of hydronephrosis
圖 3第1行左一待測樣本為一幅390*290像素的前列腺增生B超圖像，分割目標為中部團塊暗區。病灶區的大部分與背景區域的灰度差別不大，病灶內部灰度分布不勻，邊界很模糊。圖 3第1行左二的紅色閉合曲線為初始輪廓，與病灶組織的主體重合較好，但未把大部分病灶衍生的陰影組織包括進來，曲線的左右下方分別存在人工糾正的與背景斷裂的點。以手工分割結果為標準對各個模型的分割結果進行評判：在低對比度和模糊邊界的情況下，CV模型分割結果的輪廓線已經深入目標內部，并有繼續向內部區域延伸的趨勢，把一部分背景誤作為目標，直到圖像外框才停止下來，而病灶本身未完整地提取出來；RSF模型分割結果基本提取了病灶的主要部分，存在少量的雙向(背景-目標，目標-背景)誤劃分；SLGS模型分割結果提取的病灶形狀較為接近手工分割，但邊界線存在較多冗余，并且由于把大量背景劃為病灶，導致醫生所指定的病灶未形成視覺上的閉合曲線；LGF模型分割結果存在大量細碎輪廓線，視覺干擾很大；而本文方法的分割結果在形狀和面積上都很接近手工分割結果，能夠把與目標局部灰度分布概率密度接近的組織正確地劃為目標，并有效地隔離了部分偽影，但未能將增生部位完全分割出來。
圖3
前列腺增生B超圖像分割結果
Figure3.
Segmentation of BU image of prostatic hyperplasia
圖 4第1行左一待測樣本為一幅361*252像素的膽囊B超圖像，分割目標為中部橫向暗色團塊區域，炎性病灶內部區域灰度分布自上而下由稀疏到稠密，上部與背景的灰度分布呈過渡狀漸變融合，病灶邊界很模糊，左右兩側均存在面積大小不一的若干片區分布偽影，病灶內上部存在片狀噪聲。圖 4第1行左二的紅色閉合曲線為初始輪廓，與病灶組織的主體重合較好，右側有人工強制截斷點，與背景偽影相分離。以手工分割結果為標準對各個模型的分割結果進行評判：CV模型分割結果雖能大致分割出囊腫部位，但將部分病灶組織錯分為背景，并且把右側無關區域也劃歸為病灶，即使調整了參數輪廓線也比較細碎；RSF模型分割結果對病灶主體部分邊界的勾勒比較精細，但未隔離病灶與右側的偽影；SLGS模型分割結果提取的病灶形狀比RSF模型更為準確，但同樣把右側無關區域和左側部分背景也劃歸為病灶；LGF模型分割結果對病灶部分的提取僅勉強勾勒出輪廓，并產生了大量冗余細碎輪廓線；而本文方法分割結果將炎性膽囊較為準確地提取出來并形成了閉合曲線，結果最符合專業醫師的要求。
圖4
膽囊B超圖像分割結果
Figure4.
Segmentation of BU image of cholecyst
由以上實驗結果可知，CV模型借助全局信息提高了抗噪能力，但未考慮局部區域信息，水平集曲線僅根據區域均值而粗略逼近目標邊界，往往到達最大迭代步數時只能停止在邊界附近，無法準確分割灰度不勻、偽影干擾和目標背景差別小的目標，且需要對所有圖像點進行迭代計算，耗時較多；RSF模型主要依賴局部邊界梯度信息，對于邊界明顯的提取較為準確，但對模糊邊界的判斷常出現錯誤；SLGS模型綜合了CV模型和RSF模型的長處，既能夠較好地提取弱邊界，又有一定的抗噪性，并且分割速度很快，但由于采用了全局分割策略，所以無法剔除背景中與病灶相似的成分；LGF模型能夠較準確地表現局部的不規則變化，但高斯精度由核函數的半徑確定，在上述實驗中由于半徑取值較小，所以產生了大量冗余輪廓碎片；本文方法采用小窗口鄰域計算輪廓線上像素點的均值和方差，以取得其對目標和背景兩個區域的隸屬度因子，較為準確地判斷了弱邊界上點的歸屬，并且由于把計算范圍限制在初始輪廓的小鄰域，忽略了其余部分，所以計算量很小，運行時間比上述多數模型有明顯減少，同時避免了冗余輪廓的產生。
4.2 實驗結果的定量分析
下面進行分割精度的比較。本實驗以專業醫師的手工分割結果作為金標準，來評判各個模型的自動分割結果。對于分割結果存在多個連續區域的，取最大的區域或醫師手工指定的區域，從而較為客觀地驗證各模型分割的精度。醫學診斷試驗的標準化計算方法采用4個經典指標：TP(True Positives)為真陽性(實際為陽性而判斷為陽性)，FP(False Positives)為假陽性(實際為陰性而判斷為陽性)，FN(False Negatives)為假陰性(實際為陽性而判斷為陰性)，TN(True Negatives)為真陰性(實際為陰性而判斷為陰性)。在此基礎上定義本文實驗的分割精度指標[12]： TP：如果模型自動分割結果曲線內部(包括邊界)像素點在手工分割曲線內部(包括邊界)，則計數到TP中；TN：如果模型自動分割結果曲線外部像素點在手工分割曲線外部，則計數到TN中；FP：如果模型自動分割結果曲線內部(包括邊界)像素點不在手工分割曲線內部(包括邊界)，則計數到FP中；FN：如果模型自動分割結果曲線外部像素點不在手工分割曲線外部，則計數到FN中。
測準度Accuracy=，測準度又稱測準率，反映按該篩檢實驗的標準值被正確判定陽性和陰性者占總樣本量的比值。
靈敏度Sensitivity=，靈敏度又稱真陽性率，即實際有病而按該篩檢實驗的標準被正確判斷為有病的百分比，它反映篩檢實驗發現患者的能力。
特異度Specificity=，特異度又稱真陰性率，即實際無病按該診斷標準被正確地判斷為無病的百分比，它反映篩檢實驗確定非患者的能力。
表 1為5種模型處理圖 2~圖 4的迭代次數、耗時和分割精度指標測量值。
表1
5種水平集模型的迭代次數、耗時和精度比較
Table1.
Comparison of 5 level set models on iteration times,time consumption and accuracy
從南方醫院臨床醫學影像圖像庫中隨機抽取20幅B超圖像，用5種模型分別進行處理，以獲得處理精度的統計結果。由于各幅圖像大小不盡相同、內容不統一、質量不一致，所以不進行迭代次數和耗時的比較，只進行處理精度的統計比較，測量結果數據見表 2。
表2
5種水平集模型平均分割精度與是否需要后續處理
Table2.
Comparison of 5 level set models on accuracy indexes and post-processing necessity
從以上樣本統計數據可以得出結論：在分割速度方面，LGF模型迭代次數最少，耗時最少，本文模型迭代次數多于SLGS模型，少于CV模型和RSF模型，但耗時少于SLGS模型、CV模型和RSF模型。說明本文的局部計算方案具有一定的減少迭代次數、減少計算量的作用。
在測準度方面，本文模型優于其它模型，表明本文模型的陽性和陰性判斷的總體準確程度最高，CV模型、RSF模型和SLGS模型的準確率起伏較大，說明受圖像內容或質量的影響較大，而LGF模型的準確率最低；在靈敏度方面，LGF模型表現最好，本文模型次之，其他模型略遜之，表明LGF模型和本文模型較為適合發現病變組織；在特異度方面，本文模型最好，SLGS模型略次之，但起伏比本模型大，說明受圖像內容或質量的影響較大，而其他模型表現則較為遜色，表明本文模型和SLGS模型在篩檢確定非患者的能力方面較強。
5 總結
本文提出一種面向B超圖像病灶分割的水平集方法，利用隸屬度判斷像素點區域歸屬，根據輪廓線像素點對于輪廓線內部和外部的能量差值來決定移動方向，達到能量函數的最小值。通過直接計算能量函數來確定水平集的更新，不采用傳統的解偏微分方程方法，不需要掃描整個圖像，僅對輪廓曲線內外鄰域點實行轉換判斷處理，減小了計算量，較好地解決了灰度分布不勻以及低對比度B超圖像目標模糊邊界的確定問題。在與幾種優良模型的檢測比較中，本文模型在平均測準度、平均靈敏度和平均特異度方面都表現較好，并且起伏很小，不容易受圖像內容和質量的影響。
本文方法的局限性是：一、要求初始輪廓較為接近目標邊界，否則需要增大鄰域窗口的大小，從而降低速度，同時分割精度也受到很大影響。這一點決定了本文模型對手工初始化的依賴性。實際上，窗口大小變為整幅圖像時則退化為CV模型，即本文模型是CV模型的局域特例。二、本文的隸屬度判別模型針對B超灰度圖像，只利用了灰度分布這一項描述算子。根據人眼視覺注意機制的研究結果，對人眼影響較大的低層視覺顯著性特征有亮度對比度、目標區域大小、目標位置、目標形狀等，將來可根據不同內容的圖像和不同需求的分割任務進行其他描述算子的提取和設計，以實現面向具體應用需求的更好分割。
引言
超聲診斷是將超聲檢測技術應用于人體，通過測量了解生理或組織結構的數據和形態，從而發現疾病的一種診斷方法，主要用于腫物、畸形、結石及其他能引致局部結構發生明顯形態改變的疾病，具有無創、方便、廉價、實時等優點。其中應用最廣的B型超聲(B-type ultrasound,BU)檢查對液體與實質性組織具有顯著的圖像差別，通過黑白之間的不同灰度來反映不同組織回聲的強弱，對囊性液體與實質性病變組織具有良好的鑒別能力[1]。提取病灶區域是醫學診斷和制定治療計劃的先行步驟。醫學圖像分割根據區域間的相似或不同，把圖像分割成若干區域，提取感興趣組織或部位。由于超聲成像原理、被測生物體組織對聲阻抗分布不勻及在檢測過程中可能發生的運動等因素影響，B超圖像經常存在灰度不均勻、對比度較低、目標邊界模糊、偽影和噪聲干擾等現象，給自動分割帶來了較大的困難。
1 超聲分割方法的探索
1.1 分割方法概述
針對這些問題，近年來研究人員做了大量工作，所使用的方法主要有區域增長法、動態閾值分割法、模糊聚類法、分水嶺法、隱馬爾可夫模型(hidden Markov models，HMM)融合法、最大估值算法(estimation maximization，EM)和高斯混合分布模型、小波變換法等，提出了針對一些具體問題的解決方案。基于變分的水平集方法是近年來研究頗為活躍的一個分支，它將圖像分割問題表達為能量函數的最小化問題，由變分原理將其轉化為偏微分方程的求解[2]。相比于傳統的區域分割方法，變分方法可以在能量函數中加入與圖像內容有關的約束條件等綜合信息，獲得更加自然準確的分割效果。
1.2 水平集方法的應用評述
水平集方法的主要思想是將移動的界面作為零水平集嵌入高一維的水平集函數中，由封閉超曲面零水平集的演化方程得到水平集函數的演化方程，通過水平集函數曲面的演化來隱含地求解輪廓曲線，結果得到目標輪廓閉合曲線。Chan-Vese(CV)模型[3]是區域模型的代表，設圖像I(x,y)的定義域為Ω，設圖像I(x,y)被輪廓曲線C劃分為目標(曲線C的內部)和背景(曲線C的外部)兩個同質區域，平均灰度分別為Cin 和Cout，設計基于整個圖像的區域能量函數，通過面積和長度約束能量項保持曲線的連續性和平滑性，能量泛函為：
| ${E_{CV}}\left( {{c_{{\rm{in}}}},{c_{{\rm{out}}}},C} \right) = {E_{{\rm{global}}}} + {E_{{\rm{length}}}} + {E_{{\rm{area}}}} = {\lambda _1}\int\limits_{{\Omega _{{\rm{in}}}}} {{{\left| {I\left( {x,y} \right) - {c_{{\rm{in}}}}} \right|}^{\rm{2}}}{\rm{d}}x{\rm{d}}y} + {\lambda _2}\int\limits_{{\Omega _{{\rm{in}}}}} {{{\left| {I\left( {x,y} \right) - {c_{{\rm{in}}}}} \right|}^{\rm{2}}}{\rm{d}}x{\rm{d}}y} + \mu \left| C \right| + \gamma A$ |
CV模型利用圖像的全局灰度信息而非邊緣信息，能夠有效分割出模糊或部分缺失的邊緣；當兩種區域的灰度差別較大時具有較強的抗噪能力；對初始輪廓線位置不敏感；經嚴格證明能夠收斂到全局最優。但是，它的前提條件是圖像灰度分布均勻、前景背景有明顯差異，而B超圖像普遍對比度不高，分割結果不夠準確。另外，μ的取值為經驗式的，取值小時能較準確刻畫病灶不規則邊界，但輪廓線常產生大量零碎冗余，取值大時輪廓線形狀光滑卻容易偏離目標的不規則邊界。
近年來針對水平集在醫學中的運用出現了一些較好的改進方案。如文獻[4]采用了圖像平滑和圖像分割同時進行的策略，可以自動分割出多種類型的區域，能自動處理拓撲變化，但由于平滑和分割同時進行導致計算量較大，僅適用于圖像內容比較簡單的情況。文獻[5]設計了一種區域可擴展自適應能量函數(region-scalable fitting energy，RSF)以構造水平集模型，對局部區域內輪廓線的擬合能量函數捕獲圖像的局部區域信息，能夠較好地處理灰度不勻的情況，但運算代價比CV模型顯著增加，使用的核函數的方差為定值，沒有區別考慮不同區域的強度變化，并且容易陷入局部極值解。文獻[6]提出可選擇的局部或全局能量函數(selective local or global segmentation，SLGS)，運用可選擇的二值高斯過濾糾正分割方案，在弱邊界處逼停演化的曲線，具有較好的抗噪性，并且能夠較合理地提取弱邊界，但采用局部區域均值描述區域灰度仍存在一定誤差，這種情況在低對比度B超圖像中經常出現。文獻[7]提出基于局部高斯分布的自適應能量函數(local Gaussian distribution fitting energy，LGF)，將高斯核函數作為權值來刻畫局部區域灰度變化情況，能夠較準確地表現局部的不規則變化，但精度由核函數的半徑確定，不可改變。文獻[8]采用局部統計信息引導CV模型演化，思路與文獻[7]異曲同工，其優點是采用了快速水平集方法，建立存儲曲線內外鄰點的鏈表，根據曲線上點的鄰域的適應能量場對最佳逼近目標邊界位置的收斂條件進行定義，實現演化曲線只需更新這兩個表，計算代價減小了，但維系鏈表比使用數組代碼復雜度增大，精度仍然受到局部區域取值的限制。
本文針對B超圖像特點，利用區域水平集模型的優點，提出一種新的水平集方法，定義反映像素點對目標/背景兩個區域隸屬度的因子，通過概率分布估計模型計算和比較像素點的隸屬度，以此為依據對像素點進行區域劃分，以提高水平集曲線定位B超圖像模糊邊界的能力；同時本文將B超圖像目標分割看作對感興趣目標區域的局部分割，將水平集的計算求解約束到局部分割范圍，減少計算量，提高分割速度，實現病灶等特定目標區域的高效提取。
2 隸屬度區域水平集模型算法
2.1 本文高斯混合模型和水平集模型結合的思路
高斯分布(Gaussian distribution)描述了一種圍繞某個單值聚集分布的隨機變量，是統計學中最廣泛應用的一類分布。中心極限定理表明，采樣的均值近似服從高斯分布，高斯分布的信息熵在所有的已知均值及方差的連續分布中最大，這使得它成為一種均值以及方差已知的分布的自然選擇。高斯分布是一個單極值分布，不能對多極值復雜分布提供一個較好的近似，為此提出了高斯混合模型(Gaussian mixture model，GMM)，把數據看作是從數個高斯分布中生成出來的，將一個事物分解為若干的基于高斯概率密度函數(正態分布曲線)形成的模型，以得到更精確的量化模型。GMM的概率密度函數就是幾個高斯概率密度函數的加權和。圖像灰度直方圖反映的是圖像中某個灰度值出現的頻次，可以認為是圖像灰度概率密度的估計。如果圖像所包含的目標區域和背景區域在灰度上有明顯差異，那么該圖像的灰度直方圖呈現雙峰形狀，分別對應目標和背景的中心灰度。而復雜的圖像如醫學圖像、遙感圖像等，直方圖一般是多峰的。通過將直方圖的多峰特性看作是多個高斯分布的疊加，可以解決圖像的分割問題。
本文將GMM和水平集模型結合起來，定義待處理B超圖像為I(x)=a(x)r(x)，其中r(x)為理想圖像(灰度均勻分布、無噪、無偽影)，I(x)為實際待處理圖像，a(x)為偏差因子。假設理想圖像r(x)中像素點的強度分布滿足GMM：
| $\pi \left( x \right) = \sum\limits_{i = 1}^M {{\omega _i}g\left( {x;{m_i},{\sigma _i}} \right)} $ |
其中M為分割結果中的類別(區域)數目，本文取M=2(目標和背景)，ωi為先驗概率，滿足，g(x;mi,σi)為均值為mi、協方差為σi的Gaussian分布：g(x;mi,Σi)=。
式(3)表示的概率模型似然函數為
| $L\left( {u\left( x \right)} \right) = \prod\limits_{i = 1}^M {\prod\limits_{{\Omega _i}} {{\omega _i}{p_i}\left( {u\left( x \right)} \right)} } $ |
最大似然法(maximum likelihood,ML)是一種點估計法，當從模型總體隨機抽取n組樣本觀測值后，最合理的參數估計量應該使得該n組樣本觀測值的概率最大。根據最大似然估計法，求總體參數u(x)的極大似然估計值的問題就是求似然函數L(u(x))的最大值(或-L(u(x))的最小值)問題，可通過解=0來解決。因為logL是L的增函數，兩者在u(x)的同一值處取得最大值，就把問題轉化為logL的極值求解。
2.2 本文隸屬度區域水平集方法
設C={Cin,Cout}是圖像I關于輪廓線C的兩個劃分，Cin代表輪廓線內部即目標區域，Cout代表輪廓線外部即背景區域，根據文獻[9]，設p(Cin)和p(Cout)分別為目標區域和背景區域關于像素點x的先驗概率，則有p(Cin)=p(Cout)。像素點x屬于目標或背景的概率分別為
| $\begin{array}{l} p\left( {C\left. {_{{\rm{in}}}} \right|x} \right) = \frac{{p\left( {{C_{{\rm{in}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right)}}{{p\left( {{C_{{\rm{in}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right) + p\left( {{C_{{\rm{out}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right)}}\\ p\left( {C\left. {_{{\rm{out}}}} \right|x} \right) = \frac{{p\left( {{C_{{\rm{out}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right)}}{{p\left( {{C_{{\rm{in}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right) + p\left( {{C_{{\rm{out}}}}} \right)p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right)}} \end{array}$ |
其中，u(x)為I(x)中像素點x的灰度值，定義關于每個像素點x的區域隸屬度因子p(u(x)｜Cin)和p(u(x)｜Cout)，分別表示在目標類與背景類條件下出現u(x)的概率，它們服從互相獨立的高斯分布。由p(Cin)p(u(x)｜Cin)+p(Cout)p(u(x)｜Cout)為常數，p(Cin)=p(Cout)=0.5，將式(5)簡化為：
| $\begin{array}{l} p\left( {C\left. {_{{\rm{in}}}} \right|x} \right) = {\lambda _1}p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{in}}}}} \right) = \frac{\lambda }{{{\sigma _1}}}\exp \left\{ { - \frac{{{{\left( {x - {m_1}} \right)}^2}}}{{2\sigma _1^2}}} \right\}\\ p\left( {C\left. {_{{\rm{out}}}} \right|x} \right) = {\lambda _2}p\left( {u\left. {\left( x \right)} \right|{C_{{\rm{out}}}}} \right) = \frac{\lambda }{{{\sigma _2}}}\exp \left\{ { - \frac{{{{\left( {x - {m_2}} \right)}^2}}}{{2\sigma _2^2}}} \right\} \end{array}$ |
其中。
根據區域隸屬度因子進行像素點的區域歸屬判定，當p(Cin｜x)>p(Cout｜x)時，像素點x劃歸為目標類；p(Cin｜x)<p(Cout｜x)時，像素點x劃歸為背景類；p(Cin｜x)=p(Cout｜x)時，則像素點x為目標類和背景類的分界點。將目標和背景區域的符號取相反(內負外正)以實現求差運算，定義能量函數為：
| $\begin{array}{l} E\left( C \right) = - \lambda \int\limits_\Omega {p\left( {C\left. {_{{\rm{in}}}} \right|x} \right){\rm{d}}{R_{{\rm{in}}}}} + \lambda \int\limits_\Omega {p\left( {C\left. {_{{\rm{out}}}} \right|x} \right){\rm{d}}{R_{{\rm{out}}}}} + \mu \left| C \right| + \gamma A = - \lambda \int\limits_\Omega {\left[ {\log \frac{{{\sigma _1}}}{{{\sigma _2}}} - \frac{{{{\left( {u\left( x \right) - {m_1}} \right)}^2}}}{{\sigma _1^2}} + \frac{{{{\left( {u\left( x \right) - {m_2}} \right)}^2}}}{{\sigma _2^2}}} \right]} {\rm{d}}x + \\ \mu \left| C \right| + \gamma A \end{array}$ |
顯然，當p(Cin｜x)=p(Cout｜x)時式(6)達到最小值，水平集停止在目標和背景的邊界上，得到分割輪廓線。使用高斯分布計算概率分布值需要先固定分割區域，求出各區域內所有像素灰度的均值mi和方差σi作為概率分布參數。在水平集曲線演化過程中分割區域是隨輪廓線時刻變化的，每迭代一次就要根據當前的分割情況采樣所有點，計算不同區域的概率分布參數計算量很大。考慮到B超圖像中不同的組織具有不同的灰度值，每個區域內部也存在強度不均勻性，以小鄰域信息代替全局信息能夠更準確地刻畫圖像局部特征。在每個像素點x的小鄰域內計算圖像的統計信息，以更準確地反映出這種變化的強度分布。引入以像素點x為中心，大小為k×k的窗口函數wk(x)，把窗口wk(x)中曲線C的內部區域和外部區域的灰度均值分別記為minw(x)和moutw(x)，方差分別記為σinw2(x)和σoutw2(x)，內外兩個區域的像素點個數分別為ninw(x)和noutw(x)，灰度均值和方差的計算式為：
| $\begin{array}{l} {m_{{\rm{inw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{in}}}}\left( x \right)} {I\left( x \right)} }}{{{n_{{\rm{inw}}}}}},{m_{{\rm{outw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{out}}}}\left( x \right)} {I\left( x \right)} }}{{{n_{{\rm{outw}}}}}},\\ {\sigma _{{\rm{inw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{in}}}}\left( x \right)} {{{\left( {I\left( x \right) - {m_{{\rm{inw}}}}\left( x \right)} \right)}^2}} }}{{{n_{{\rm{inw}}}}}},\\ {\sigma _{{\rm{outw}}}}\left( x \right) = \frac{{\sum\limits_{{w_{k,{\rm{out}}}}\left( x \right)} {{{\left( {I\left( x \right) - {m_{{\rm{outw}}}}\left( x \right)} \right)}^2}} }}{{{n_{{\rm{outw}}}}}} \end{array}$ |
引入水平集表示法，得到隸屬度區域水平集能量演化方程：
| $\frac{{\partial \Phi }}{{\partial t}}\delta \left( \Phi \right)\left[ {\log \frac{{{\sigma _{{\rm{inw}}}}}}{{{\sigma _{{\rm{outw}}}}}} - \frac{{{{\left( {u\left( x \right) - H\left( {\Phi \left( {u\left( x \right)} \right)} \right) * {m_{inw}}} \right)}^2}}}{{\sigma _{{\rm{inw}}}^2}} + \frac{{{{\left( {u\left( x \right) - \left( {1 - H\left( {\Phi \left( {u\left( x \right)} \right)} \right)} \right) * {m_{{\rm{outw}}}}} \right)}^2}}}{{\sigma _{{\rm{outw}}}^2}} + \mu {\rm{div}}\left( {\frac{{\nabla \Phi }}{{\left| {\nabla \Phi } \right|}}} \right) + \gamma } \right]$ |
其中Hε(x)為Heaviside函數，記做Hε(x)=，δε(ф)=H′ε(x)=。這樣就把圖像目標分割問題轉化為極小化能量函數(8)的計算。
3 計算方法的具體實現
3.1 初始輪廓的提取
大多數區域水平集模型是全局分割模型，將圖像中的所有目標分割出來，計算量很大，對于目標針對性強的B超圖像無此必要。本文水平集模型對初始輪廓的設置比較敏感，初始輪廓的設置得當，有助于提高分割準確度和效率。超聲病灶在不同圖像中位置并不固定，如果將初始輪廓統一設置在固定位置、統一形狀(如橢圓或矩形)，會導致初始輪廓與目標存在較大偏離，造成計算時間過長，甚至產生錯誤的分割結果。
經過比較測試，本文采用計算簡單的最大類間方差法(Otsu算法)自動確定圖像分割的初始條件。最佳閾值法按圖像的灰度特性將圖像分成背景和前景兩部分，背景和前景之間的類間方差越大，說明構成圖像的兩部分的差別越大，當部分前景錯分為背景或部分背景錯分為前景時都會導致兩部分差別變小，設灰度圖像灰度級是L，灰度范圍為[0,L-1]，最佳閾值計算公式為：
| $t = \max \left( {{w_b}\left( t \right) * {{\left( {{u_b}\left( t \right) - \bar u} \right)}^2} + {w_f}\left( t \right) * {{\left( {{u_f}\left( t \right) - \bar u} \right)}^2}} \right)$ |
其中，wb為背景比例，ub為背景均值，wf為前景比例，uf為前景均值，
Otsu二值化結果中常存在著許多碎片，需要進行后續處理，本文將數學形態學的開閉運算結合起來對二值圖像進行噪聲濾除，以開運算為主，先腐蝕后膨脹，以消除孤立噪音等小物體、斷裂背景與目標的纖細連接處、平滑較大物體的邊界；以閉運算為輔，先膨脹后腐蝕，填充目標內部細小空洞、連接鄰近物體、平滑邊界。在此基礎上運用八鄰域跟蹤算法勾勒病灶邊界，提取初始輪廓，作為水平集的起點進行演化運算。
圖 1為提取初始輪廓的過程。需要注意的是，由于部分B超圖像內容復雜，當出現多個閉合曲線時需要人工選擇目標區域；當閉合曲線與背景有錯誤粘連時，需要人工截斷錯處的連接點，設置正確的起點，重新進行初始輪廓點的搜索。
圖1
確定水平集初始輪廓的過程
(a)待處理素材；(b)Otsu二值化結果；(c)數學形態學開閉運算的濾波結果；(d)八鄰域跟蹤算法提取的閉合曲線
Figure1. Procedure of getting initial contour for level set(a) image to be processed； (b) result of Otsu binary process； (c) result of mathematical morphology filting； (d) closed curve achieved by eight-neighborhood tracking
3.2 快速高效的局部計算方案
為提高計算效率，本文水平集模型采用局部計算方案，定義基于局部近似符號距離函數，放寬對初始水平集函數的要求，不再要求初始水平集必須是符號距離函數。采用曲線上每個像素點周圍的k×k個鄰居點的聯合概率分布代替該點的概率分布，以減少由于灰度分布不均、對比度低和噪聲干擾帶來的計算誤差，比采用全局信息更能準確地定位模糊邊界。但如果直接對輪廓線上的每個像素點都做對數及求和運算，求解水平集式需要離散化迭代實現求解偏微分方程，計算量較大，曲線演化速度較慢。因此采用局部計算思路以減少運算量和運算時間。文獻[9-10]提出了改進的計算方案以及免除反復初始化水平集的做法，能夠有效加快演化速度，但某些實現技術較為復雜。本文借鑒文獻[11]，采用改進的FTC(Fast two-cycle)算法，在不用求解偏微分方程的基礎上近似實現水平集函數的進化，更新局部近似符號距離函數的定義為
| ${\Phi _0}\left( {x,y} \right) = \left\{ {\begin{array}{*{20}{l}} {k,點x為外部點}\\ {1,點x為輪廓線外領域點}\\ {0,點x在輪廓線上}\\ { - 1,x為輪廓線內領域點}\\ { - k,點x為內部點} \end{array}} \right.$ |
水平集迭代過程是：
(1)如果輪廓線向外運動一個像素，則判斷Lout中的每個點x的能量變化ΔELout(新能量值－原能量值)的符號情況：
①如果ΔELout≥0則不處理；
②如果ΔELout<0，則輪廓線在該點向外擴展一個像素，并且掃描Lin，如果存在x的鄰接點y，則修改ф(y)為-k，即當曲線向外擴展時將原來的內鄰域點修改為內部點。
(2)如果輪廓線向內運動一個像素，則判斷Lin中的每個點x的能量變化 ΔELin(新能量值-原能量值)的符號情況：
①如果ΔELin≥0則不處理；
②如果ΔELin<0，則輪廓線在該點向內擴展一個像素，并且掃描Lout，如果存在x的鄰接點y，則修改ф(y)為k，即當曲線向內擴展時將原來的外鄰域點修改為外部點。
一般情況下，k值取得較大，則分割較為精確，但考慮到計算代價，k值一般取為3或5。每次迭代結束后所得目標輪廓作為下一次迭代的初始輪廓，對圖像進行連續分割，當能量函數的變化小于一定的閾值時，迭代終止，得到分割結果。只需對輪廓線上的點而不是整幅圖像的點進行計算，通過計算能量函數的增減變化來確定輪廓線的擴展或收縮。本文方法的局限是需要保持初始輪廓接近目標邊界，否則就必須擴大鄰域窗口半徑來保證輪廓線形變的正確性，而加大鄰域窗口則意味著運算量的增加。
4 結果
分割實驗待測圖像來自南方醫科大學附屬醫院影像庫，通過對相同B超圖像樣本進行分割來檢驗本文模型的效果，并與幾種性能優良的水平集模型進行比較。編程平臺為Matlab7，計算機配置為Intel? CoreTM(2) Duo CPU，3.0 GHz，2 GB RAM，Windows XP 操作系統。所有實驗采用相同的終止閾值ΔE=10-4或最大迭代次數2 000(考慮到用戶等待時間和實用性)。圖 2~圖 4分別為臨床B超待分割圖像、本文模型初始輪廓、CV模型[3]、RSF模型[5]、SLGS模型[6]、LGF模型[7]、本文模型分割結果和專業醫師手工分割結果。所有模型均按原文獻方法進行相應參數設置，其中公共參數：時間步長Δt=0.1，空間步長Δh=1，對Heaviside函數和Dirac函數進行近似時正則化參數ε=1，長度懲罰項參數μ、鄰域半徑等則根據圖像內容調整為最合適值，以盡可能達到每種模型最優的分割效果。
從分割結果的視覺效果和統計信息兩個方面來比較5種模型在處理B超圖像方面的適用性。
4.1 實驗結果的定性分析
首先看分割結果的視覺效果。圖 2第1行左一為待測樣本411*290像素腎盂積水B超截圖，分割目標為中部位置的暗色積液區域。整幅圖像存在明顯的灰度分布不勻，有多處偽影存在，左側存在大面積偽影，上部和右側存在帶狀偽影，存在散布的噪聲雜質，積液區域的上部有兩處凸起，邊界較為模糊，左側和右側與偽影的交界比較模糊。圖 2第1行左二中的紅色閉合曲線為初始輪廓，并未將病灶完全包括在內，而是與積液區域交疊。以手工分割結果為標準對各個模型的分割結果進行評判：CV模型分割結果較好地處理了灰度不勻的現象，并具有一定的抗噪聲能力；RSF模型分割結果顯示，曲線較容易越過邊界向病灶內部運動，最終收斂到噪聲處，無法準確貼近病灶邊界；SLGS模型分割結果對不規則處的處理較好，但容易越過病灶邊界深入到病灶內部，使得病灶有效提取區比實際略小，同時也把圖像中無關的部分中與病灶灰度分布相似的成分也分割出來，造成不必要的干擾；LGF模型分割結果雖能基本準確地提取病灶邊界，但產生了大量無用的分割碎片，分割結果必須經過后期處理才能用于輔助診斷；而本文方法分割結果在病灶提取方面接近CV模型，對不規則凸起和凹陷處的提取則更為合理，在刻畫模糊邊界處也更為精細準確，能夠反映病灶形態的變化，最接近手工分割結果，但提取的病灶區比實際略大。
圖2
腎盂積水B超圖像分割結果
Figure2.
Segmentation of BU image of hydronephrosis
圖 3第1行左一待測樣本為一幅390*290像素的前列腺增生B超圖像，分割目標為中部團塊暗區。病灶區的大部分與背景區域的灰度差別不大，病灶內部灰度分布不勻，邊界很模糊。圖 3第1行左二的紅色閉合曲線為初始輪廓，與病灶組織的主體重合較好，但未把大部分病灶衍生的陰影組織包括進來，曲線的左右下方分別存在人工糾正的與背景斷裂的點。以手工分割結果為標準對各個模型的分割結果進行評判：在低對比度和模糊邊界的情況下，CV模型分割結果的輪廓線已經深入目標內部，并有繼續向內部區域延伸的趨勢，把一部分背景誤作為目標，直到圖像外框才停止下來，而病灶本身未完整地提取出來；RSF模型分割結果基本提取了病灶的主要部分，存在少量的雙向(背景-目標，目標-背景)誤劃分；SLGS模型分割結果提取的病灶形狀較為接近手工分割，但邊界線存在較多冗余，并且由于把大量背景劃為病灶，導致醫生所指定的病灶未形成視覺上的閉合曲線；LGF模型分割結果存在大量細碎輪廓線，視覺干擾很大；而本文方法的分割結果在形狀和面積上都很接近手工分割結果，能夠把與目標局部灰度分布概率密度接近的組織正確地劃為目標，并有效地隔離了部分偽影，但未能將增生部位完全分割出來。
圖3
前列腺增生B超圖像分割結果
Figure3.
Segmentation of BU image of prostatic hyperplasia
圖 4第1行左一待測樣本為一幅361*252像素的膽囊B超圖像，分割目標為中部橫向暗色團塊區域，炎性病灶內部區域灰度分布自上而下由稀疏到稠密，上部與背景的灰度分布呈過渡狀漸變融合，病灶邊界很模糊，左右兩側均存在面積大小不一的若干片區分布偽影，病灶內上部存在片狀噪聲。圖 4第1行左二的紅色閉合曲線為初始輪廓，與病灶組織的主體重合較好，右側有人工強制截斷點，與背景偽影相分離。以手工分割結果為標準對各個模型的分割結果進行評判：CV模型分割結果雖能大致分割出囊腫部位，但將部分病灶組織錯分為背景，并且把右側無關區域也劃歸為病灶，即使調整了參數輪廓線也比較細碎；RSF模型分割結果對病灶主體部分邊界的勾勒比較精細，但未隔離病灶與右側的偽影；SLGS模型分割結果提取的病灶形狀比RSF模型更為準確，但同樣把右側無關區域和左側部分背景也劃歸為病灶；LGF模型分割結果對病灶部分的提取僅勉強勾勒出輪廓，并產生了大量冗余細碎輪廓線；而本文方法分割結果將炎性膽囊較為準確地提取出來并形成了閉合曲線，結果最符合專業醫師的要求。
圖4
膽囊B超圖像分割結果
Figure4.
Segmentation of BU image of cholecyst
由以上實驗結果可知，CV模型借助全局信息提高了抗噪能力，但未考慮局部區域信息，水平集曲線僅根據區域均值而粗略逼近目標邊界，往往到達最大迭代步數時只能停止在邊界附近，無法準確分割灰度不勻、偽影干擾和目標背景差別小的目標，且需要對所有圖像點進行迭代計算，耗時較多；RSF模型主要依賴局部邊界梯度信息，對于邊界明顯的提取較為準確，但對模糊邊界的判斷常出現錯誤；SLGS模型綜合了CV模型和RSF模型的長處，既能夠較好地提取弱邊界，又有一定的抗噪性，并且分割速度很快，但由于采用了全局分割策略，所以無法剔除背景中與病灶相似的成分；LGF模型能夠較準確地表現局部的不規則變化，但高斯精度由核函數的半徑確定，在上述實驗中由于半徑取值較小，所以產生了大量冗余輪廓碎片；本文方法采用小窗口鄰域計算輪廓線上像素點的均值和方差，以取得其對目標和背景兩個區域的隸屬度因子，較為準確地判斷了弱邊界上點的歸屬，并且由于把計算范圍限制在初始輪廓的小鄰域，忽略了其余部分，所以計算量很小，運行時間比上述多數模型有明顯減少，同時避免了冗余輪廓的產生。
4.2 實驗結果的定量分析
下面進行分割精度的比較。本實驗以專業醫師的手工分割結果作為金標準，來評判各個模型的自動分割結果。對于分割結果存在多個連續區域的，取最大的區域或醫師手工指定的區域，從而較為客觀地驗證各模型分割的精度。醫學診斷試驗的標準化計算方法采用4個經典指標：TP(True Positives)為真陽性(實際為陽性而判斷為陽性)，FP(False Positives)為假陽性(實際為陰性而判斷為陽性)，FN(False Negatives)為假陰性(實際為陽性而判斷為陰性)，TN(True Negatives)為真陰性(實際為陰性而判斷為陰性)。在此基礎上定義本文實驗的分割精度指標[12]： TP：如果模型自動分割結果曲線內部(包括邊界)像素點在手工分割曲線內部(包括邊界)，則計數到TP中；TN：如果模型自動分割結果曲線外部像素點在手工分割曲線外部，則計數到TN中；FP：如果模型自動分割結果曲線內部(包括邊界)像素點不在手工分割曲線內部(包括邊界)，則計數到FP中；FN：如果模型自動分割結果曲線外部像素點不在手工分割曲線外部，則計數到FN中。
測準度Accuracy=，測準度又稱測準率，反映按該篩檢實驗的標準值被正確判定陽性和陰性者占總樣本量的比值。
靈敏度Sensitivity=，靈敏度又稱真陽性率，即實際有病而按該篩檢實驗的標準被正確判斷為有病的百分比，它反映篩檢實驗發現患者的能力。
特異度Specificity=，特異度又稱真陰性率，即實際無病按該診斷標準被正確地判斷為無病的百分比，它反映篩檢實驗確定非患者的能力。
表 1為5種模型處理圖 2~圖 4的迭代次數、耗時和分割精度指標測量值。
表1
5種水平集模型的迭代次數、耗時和精度比較
Table1.
Comparison of 5 level set models on iteration times,time consumption and accuracy
從南方醫院臨床醫學影像圖像庫中隨機抽取20幅B超圖像，用5種模型分別進行處理，以獲得處理精度的統計結果。由于各幅圖像大小不盡相同、內容不統一、質量不一致，所以不進行迭代次數和耗時的比較，只進行處理精度的統計比較，測量結果數據見表 2。
表2
5種水平集模型平均分割精度與是否需要后續處理
Table2.
Comparison of 5 level set models on accuracy indexes and post-processing necessity
從以上樣本統計數據可以得出結論：在分割速度方面，LGF模型迭代次數最少，耗時最少，本文模型迭代次數多于SLGS模型，少于CV模型和RSF模型，但耗時少于SLGS模型、CV模型和RSF模型。說明本文的局部計算方案具有一定的減少迭代次數、減少計算量的作用。
在測準度方面，本文模型優于其它模型，表明本文模型的陽性和陰性判斷的總體準確程度最高，CV模型、RSF模型和SLGS模型的準確率起伏較大，說明受圖像內容或質量的影響較大，而LGF模型的準確率最低；在靈敏度方面，LGF模型表現最好，本文模型次之，其他模型略遜之，表明LGF模型和本文模型較為適合發現病變組織；在特異度方面，本文模型最好，SLGS模型略次之，但起伏比本模型大，說明受圖像內容或質量的影響較大，而其他模型表現則較為遜色，表明本文模型和SLGS模型在篩檢確定非患者的能力方面較強。
5 總結
本文提出一種面向B超圖像病灶分割的水平集方法，利用隸屬度判斷像素點區域歸屬，根據輪廓線像素點對于輪廓線內部和外部的能量差值來決定移動方向，達到能量函數的最小值。通過直接計算能量函數來確定水平集的更新，不采用傳統的解偏微分方程方法，不需要掃描整個圖像，僅對輪廓曲線內外鄰域點實行轉換判斷處理，減小了計算量，較好地解決了灰度分布不勻以及低對比度B超圖像目標模糊邊界的確定問題。在與幾種優良模型的檢測比較中，本文模型在平均測準度、平均靈敏度和平均特異度方面都表現較好，并且起伏很小，不容易受圖像內容和質量的影響。
本文方法的局限性是：一、要求初始輪廓較為接近目標邊界，否則需要增大鄰域窗口的大小，從而降低速度，同時分割精度也受到很大影響。這一點決定了本文模型對手工初始化的依賴性。實際上，窗口大小變為整幅圖像時則退化為CV模型，即本文模型是CV模型的局域特例。二、本文的隸屬度判別模型針對B超灰度圖像，只利用了灰度分布這一項描述算子。根據人眼視覺注意機制的研究結果，對人眼影響較大的低層視覺顯著性特征有亮度對比度、目標區域大小、目標位置、目標形狀等，將來可根據不同內容的圖像和不同需求的分割任務進行其他描述算子的提取和設計，以實現面向具體應用需求的更好分割。
