LATE水平集圖像分割模型的矩形窄帶法①

曾笑云，楊晟院，潘園園，劉 洋，左國才

1(湘潭大學 信息工程學院，湘潭 411105)

2(湖南軟件職業學院 軟件與信息工程學院，湘潭 411100)

水平集最先由Osher[1]提出，是一種將低維問題嵌入高維問題求解的方法.水平集廣泛應用于圖像分割[2，3]，它可獲得亞像素精度的封閉輪廓和區域.以2D 灰度圖像二相水平集為例，水平集取值的正負將圖像域劃分為目標和背景兩個區域.水平集的零等高線被視為目標和背景的分界線，即二維的活動輪廓曲線.活動輪廓隨著水平集的演化不斷逼近真實輪廓，最終完成分割.

早期的水平集模型，如MS 模型[4]，CV 模型[5]，SLGS模型[6]，對灰度不均勻的圖像分割效果不理想.為了有效分割灰度不均勻圖像，水平集模型朝著越來越復雜的方向發展，如RSF 模型[7]，LIC 模型[8]，LSACM 模型[9]，LATE 模型[10].其中，LATE 模型利用泰勒展開式對擬合函數進行非線性逼近，極大的提高了分割灰度嚴重不均勻圖像的能力.

水平集方法雖然能獲得較好分割效果，但也提高了計算復雜度.為了提高運算效率，Chopp[11]提出窄帶法，并由Adalsteinsson 等[12]給出了詳細的實現方法.窄帶法的核心思想就是把計算區域約束到活動輪廓附近的帶狀區域，避免了對整個圖像域進行計算，以此來提高水平集方法的計算效率.

窄帶方法能夠提高水平集的效率，首先要求窄帶生成過程要盡可能快，生成窄帶新增的計算代價應小于由于窄帶減少計算面積而節約的計算代價.為了提高窄帶生成速度，目前的方法主要有快速進行法[13，14]，快速掃描法[15，16]，DTM 窄帶法[17]等.

文獻[13]認為可設置較寬的窄帶.當活動輪廓未達到窄帶邊緣時，無需更新窄帶，當活動輪廓達到窄帶邊緣，但不發生波動時，說明活動輪廓已經達到圖像邊緣，也無需更新窄帶.只有當窄帶邊界點上有活動輪廓波動變化時，才需更新窄帶.

文獻[18，19]拋棄了窄帶更新策略，在每次水平集迭代中均更新窄帶.這樣可規避繁瑣的窄帶更新條件判斷.在不損害分割質量的前提下，設置盡可能窄的窄帶，能更好地發揮窄帶縮小計算區域的作用.

文獻[20]提出窄帶壓縮數據結構.為了保留鄰域信息，分別做行方向和列方向兩次壓縮.壓縮的窄帶結構雖然可以規避對非窄帶區域的范圍判斷，但使得鄰域結構變得復雜.

雖然傳統窄帶法有效地減少了計算范圍，提高了計算效率，但傳統窄仍然存在冗余的計算區域.在水平集演化的過程中，一部分活動輪廓先到達圖像邊緣不再運動，另一部分活動輪廓還需繼續演化逼近圖像邊緣.先到達圖像邊緣的活動輪廓在后續演化過程中的計算屬于冗余計算，還未達到圖像邊緣的活動輪廓才真正需要進行演化計算.為此，我們提出活動約束策略，將窄帶的范圍進一步約束在未達到圖像邊緣的活動輪廓的區域.

約束的活動輪廓區域形狀不規則，可采用最小矩形覆蓋不規則的窄帶區域，從而構造矩形窄帶.為了保證矩形窄帶的總面積盡可能小，對窄帶區域進行了合并優化.相比傳統不規則窄帶，矩形窄帶結構更簡單，更有利于演化計算.

為此，本文提出了一種基于LATE 水平集圖像分割模型的矩形窄帶法.本文剩余部分結構安排如下：第1 節介紹水平集和傳統窄帶等相關知識；第2 節介紹本文的矩形窄帶法；第3 節為數值實驗；第4 節為結論.

1 相關知識

1.1 水平集

設 Ω ∈R2為圖像域，I:Ω →R為給定的灰度圖像.φ:Ω →R為 水平集函數[21].活動輪廓曲線C為φ 的零等高線，如圖1所示.

圖1 水平集的原理示意圖

水平集方法可概括為：首先確立合理的水平集微分方程? φ/?t并給出恰當的初始水平集.然后利用水平集微分方程不斷更新水平集，驅動活動輪廓C向圖像的真實輪廓運動，直到完成分割.

以LATE 水平集模型[10]為例，它的水平集微分方程為：

其中，μ，ν 均為實數常量，δ (·) 為 單位沖擊函數，kσ為高斯內核函數.LIMi，Ci，b的表達式分別為：

其中，y為以x為 中心的鄰域，Ri(i∈{1,2})分別表示目標區域和背景區域.

其中，M1(φ)=H(φ)，M2(φ)=1-H(φ)，H(·)為單位階躍函數.

水平集方法中，稱每次迭代的計算復雜度為水平集復雜度，記為Cls.LATE 模型中，復雜度最高的運算是卷積運算.假設圖像寬度和長度均為N，卷積模板寬度為D，其中，D遠遠小于N.則LATE 的水平集復雜度為：

由于光照磁場或成像裝置缺陷等因素[7]，可能導致獲得的圖像是灰度不均勻的.灰度不均勻給圖像分割帶來很大的挑戰.在文獻[8]中，灰度不均勻圖像被看作灰度均勻圖像與偏置光源場的混合.蝴蝶圖像與嚴重不均勻偏置場的混合如圖2所示.

圖2 灰度嚴重不均勻圖像的合成

不同水平集模型對灰度嚴重不均勻蝴蝶圖像的分割結果如圖3所示.可見，LATE 模型具有很強的分割灰度不均勻圖像的能力.但是，LATE 模型的計算代價較高，分割緩慢.

圖3 不同水平集模型針對灰度嚴重不均勻圖像的分割結果.

1.2 傳統窄帶及其復雜度分析

雖然水平集能夠獲得較好的分割結果，但它將問題提高了一個維度，增加了計算復雜度.直觀地看，活動輪廓C受過零點附近的水平集 φ值的變化影響較大，受遠離過零點的水平集 φ值影響較小.傳統窄帶法通過把整個圖像域的水平集的計算區域限制在過零點附近的窄帶上來減小計算量，如圖4所示.

圖4 窄帶的原理示意圖

假設傳統窄帶的寬度與卷積模板寬度具有相同規模，傳統窄帶長度規模與圖像寬度相當.則傳統窄帶的面積規模可以記為：Snb=O(ND).本文中傳統窄帶法泛指窄帶面積規模為Snb=O(ND)的窄帶法，如直接生成窄帶法，快速進行法[13，14]，快速掃描法[15，16]和DTM[17]窄帶法等.

窄帶水平集方法每次迭代的復雜度稱為總復雜度，記為Call.總復雜度可以分成兩部分：生成窄帶的計算復雜度稱為窄帶復雜度，記為Cnb；對應窄帶區域的水平集復雜度Snb/N2·Cls.總復雜度可表示為：

一般地，生成窄帶的復雜度Cnb不應高于水平集復雜度Cls，否則總復雜度反而上升.因此，窄帶水平集方法的總復雜度主要取決于窄帶面積的規模Snb.

當Snb固 定時，生成窄帶復雜度Cnb越小越好.直接生成窄帶法對每個像素點進行掃描，判斷其與過零點的距離是否小于半個窄帶寬度，如果是則將其加入窄帶.因此，直接生成窄帶法的生成窄帶復雜度為：=O(N2D2).DTM 窄帶法則遍歷每個過零點，在半個窄帶寬度范圍內標記窄帶.因此，DTM 的生成窄帶復雜度為：=O(ND2).

快速進行法，快速掃描法以及DTM 方法等窄帶法都只是降低了生成窄帶復雜度Cnb.由于窄帶水平集方法的總復雜度主要取決于窄帶面積的規模Snb.因此，傳統窄帶法對提高水平集分割圖像的效率有限.

2 本文的矩形窄帶法

2.1 尋找過零點

設i,j分別表示2D 圖像域上x軸和y軸上的坐標值，I(i,j) 為圖像，φ (i,j) 為水平集，C為活動輪廓，則C為φ(i,j)的零等高線.

設CRS(i,j) 為二值矩陣，CRS(i,j)=1表示像素點(i,j) 為過零點，CRS(i,j)=0 表示像素點(i,j)為非過零點.VIS(i,j) 為 訪問矩陣，用來標記某像素點(i,j)被標記為過零點的次數.初始時CRS與VIS均為零矩陣，每次更新水平集CRS重新置零.可在x軸和y軸正負4 個方向尋找過零點.對每個像素點(i,j)做以下操作：

則可找到當前所有的過零點CRS(i,j)=1，以及統計當前像素點被訪問的次數VIS.

在水平集演化的過程中，一部分活動輪廓先到達圖像邊緣不再運動，另一部分活動輪廓還需繼續演化逼近圖像邊緣.先到達圖像邊緣的活動輪廓在后續演化過程中的計算屬于冗余計算，還未達到圖像邊緣的活動輪廓才真正需要進行演化計算.如果在靜止的過零點附近生成窄帶，對活動輪廓的演化并沒有作用.可以設置延時參數delay，若某像素點被訪問的次數超過delay次，則認為它是靜止不動的，將不再運動的過零點從CRS(i,j)中剔除.稱這樣的處理方法為活動約束.對每個像素點(i,j)做活動約束，具體操作如下：

在活動約束之前，過零點長度的規模為N.假設delay的規模與窄帶寬相同，則活動約束處理之后，可認為過零點的長度規模為D.

2.2 生成矩形窄帶

傳統窄帶的生成，是在過零點基礎上通過偏移生成帶狀區域的過程.傳統方法獲得的過零點曲線是封閉的，而經過活動約束，本文方法獲得的過零點集合只是傳統過零點曲線的一部分，是多段開放的鏈，稱為過零點鏈.原則上，應該在每條過零點鏈的鄰域內生成不規則的窄帶.由于活動約束下每條過零點鏈長度較短，窄帶區域接近矩形區域，不妨用規則的矩形區域代替不規則的窄帶區域.傳統的窄帶區域是不規則的，需要單獨建立窄帶中每個點與原圖像域之間的映射，才能將窄帶的計算結果返回到圖像域.而矩形區域則只需給出矩形的兩個對角端點便可確定一個區域，且矩形域與圖像域的像素點位置對應非常簡單，這更有利于快速卷積運算以及邊界處理.

可利用掃描法獲得每條過零點鏈對應的矩形窄帶區域，步驟如下：

1)分別對VIS中每個值為1 的像素點作為掃描的起點，逐一進行掃描.將該點坐標插入鏈碼q，將該像素點的VIS值置零.

2)若在8 鄰域內掃描到VIS為1 的點，將該點作為新的掃描起點，將其坐標點插入鏈碼q，對應VIS值置零.

3)重復步驟2)，直到在8 鄰域內找不到下一個VIS為1 的點.

4)確定矩形域.取鏈碼q中最小的橫坐標和縱坐標作為矩形框的第一個端點，取最大的橫坐標和縱坐標作為矩形框的第二個端點.得到最小的覆蓋矩形域.將矩形區域向四周擴寬大D/2 個寬度，防止過零點暴露在窄帶最外層.

活動約束后，剩余過零點的長度規模為O (D).生成矩形窄帶的過程既是對活動約束后剩余過零點的掃描過程，該過程不再對半個窄帶寬度范圍內的點進行判斷.因此，本文的方法生成窄帶的計算復雜度僅為=O(D) .矩形窄帶的長寬規模均為O (D)，因此，窄帶的面積規模僅為=O(D2).

尋找過零點方法具有簡單高效的優點，但可能存在個別過零點的遺落，導致過零點曲線不連續，從而產生面積重疊的矩形區域.如圖5中的左圖所示，矩形A與矩形B產生重疊.記矩形A的面積為SA，矩形B的面積為SB.對于重疊的矩形A和矩形B，將右上和左下的區域補充形成一個大的矩形C，如圖5所示，記矩形C的面積為SC.如果s，則矩形A和矩形B不需合并成矩形C.反之，則采用矩形C替換矩形A和矩形B.這樣處理可以減少計算區域，從而提高計算效率.

圖5 矩形區域合并優化

具體操作如下：設 (sx1,sy1)，(ex1,ey1)分別表示矩形框A的兩個端點，(sx2,sy2)，(ex2,ey2)分別表示矩形框B的兩個端點，如圖5所示.

則判斷它們相交的條件為：

令x1=min(sx1,ex1)，y1=min(sy1,ey1)，x1=min(sx1,ex1)，y2=min(sy2,ey2).判斷合并后面積更小的條件為：

每個矩形其實可以看做一個已知兩個端點的二維閉區間，面積優化只需對閉區間的端點進行操作，計算代價遠遠小于尋找的遺落過零點的代價.

在每次水平集迭代中，我們生成窄帶的算法總結如算法1 所示.

由于圖像也是矩形數組，矩形框的窄帶結果和原圖像能保持一致，直接將每個矩形框代入原水平集函數計算即可.返回結果時只需做坐標偏移處理，便可實現窄帶到原圖像的映射，從而在本次迭代中實現水平集的更新.可見，本文的方法與水平集方法非常容易實現對接.

2.3 與傳統窄帶法的比較

傳統窄帶示意圖如圖6(a)所示.傳統窄帶存在很大部分的無效計算區域，實際的活動窄帶區域只占總窄帶的一小部分.以DTM 方法為例，DTM 方法生成窄帶的復雜度為=O(ND2)，窄帶面積規模為=O(ND).根據式(5)，DTM 方法的總復雜度為：

本文的矩形窄帶示意圖如圖6(b)所示.通過活動約束，對長期禁止不動的點進行了屏蔽，大幅減少了窄帶的生成范圍.掃描生成矩形窄帶的復雜度為：=O(D)，窄帶面積規模為=O(D2).根據式(5)，矩形窄帶法的總復雜度為：

由于D遠遠小于N，所以＜可見，本文的矩形窄帶方法不僅生成窄帶復雜度低于DTM 方法，而且結合LATE 水平集的總復雜度也低于DTM 方法.其中，窄帶面積規模的降低，對進一步提高窄帶效率起到關鍵作用.

圖6 矩形窄帶和傳統窄帶示意圖

3 數值實驗

3.1 窄帶面積比較

本文的實驗結果均在MATLAB R2016a 上實現，操作系統為Win10.圖7(a)(b)分別為針對不同灰度不均勻情形的兩組實驗.實驗采用矩形窄帶法結合LATE模型[10]生成水平集，并利用水平集過零點分別生成傳統窄帶(TradNb)、僅添加活動約束的不規則窄帶(AcNb)以及矩形窄帶(RecNb).

圖7 窄帶演化過程.(a)從左至右迭代次數分別為1、44、88、132、176、220、266；(b)從左至右迭代次數分別為1、25、50、75、100、125、150.第1 行：圖像的活動輪廓；第2 行：傳統窄帶(TradNb)；第3 行：僅添加活動約束的不規則窄帶(AcNb，本文提出的過渡方案)；第4 行：活動約束矩形窄帶(RecNb，本文最終采用方案).

本文中傳統窄帶法(TradNb)泛指窄帶面積規模為Snb=O(ND)的一類窄帶法，如直接生成窄帶法，快速進行法，快速掃描法和DTMP 窄帶法等.AcNb 為本文提出的過渡方案，RecNb 為本文最終采用方案.根據前文分析，不規則窄帶矩形化的過程可能略微增加窄帶區域.因此，RecNb 相比AcNb 面積可能略微增加，但AcNb 的計算和實現更為簡單.RecNb 和AcNb 的面積規模均為Snb=O(D2)，其中D＜N.

圖7(a)為灰度不均勻程度一般的圖像，窄帶半徑為5.圖7(b) 為灰度嚴重不均勻的圖像，窄帶半徑為10.LATE 利用泰勒展開式對灰度不均勻進行調節.當灰度變化平緩時，各像素點的水平集受其鄰域的影響較小.當灰度不均勻程度較嚴重時，較遠鄰域的泰勒展開權值變大，LATE 模型能夠利用較遠鄰域信息對灰度不均勻進行修正.因此，在灰度嚴重不均勻區域，LATE水平集分割緩慢且更容易受到窄帶的影響.可知，灰度不均勻程度越高，所需要設置的最小窄帶半徑越大.

圖7中第3 行的窄帶面積明顯小于第2 行，表明活動約束能夠很大程度地減少窄帶范圍.RecNb 為AcNb 的最小矩形區域，RecNb 與AcNb 的面積規模相當，但矩形區域更方便計算機處理.圖7中第3 行與第4 行表明RecNb 與AcNb 的窄帶位置和面積相差不大.

圖7中TradNb、AcNb 和RecNb 對應的窄帶面積如表1和表2所示.可見，AcNb 的窄帶不大于TradNb的窄帶面積.理論上，RecNb 的面積略大于AcNb 的面積.但在實際掃描生成矩形窄帶的過程中，一些不必要的過零點被拋棄，RecNb 的面積也可能略小于AcNb.從表1和表2的數據來看，RecNb 與AcNb 的窄帶面積相差不大.圖8表明隨著水平集的演化，矩形窄帶的面積與傳統窄帶面積之比逐漸減少到0.可見，本文的矩形窄帶法能有效地減少窄帶面積，從而提高計算效率.

表1 圖7(a)對應窄帶面積

表2 圖7(b)對應的窄帶面積

圖8 圖7中的矩形窄帶面積與傳統窄帶面積的比值

3.2 運行效率比較

為了表述方便，本節將直接窄帶與LATE 模型結合的窄帶水平集記為DRCTLS；DTM 窄帶與LATE 模型結合的窄帶水平集記為DTMLS；矩形窄帶與LATE 模型結合的窄帶水平集記為RECLS.LATE 方法，DRCTLS 方法，DTMLS 方法以及RECLS 方法對圖像的分割結果對比如圖9所示.在對比實驗中，保證各方法的共有的參數完全相同.

當灰度不均勻程度一般時(圖9(a)-圖9(e))，各方法都能得到較好的分割結果.當灰度嚴重不均勻時，灰度不均勻區域分割緩慢，且容易受到窄帶范圍的影響.在相同的窄帶寬度下，DRCTLS 方法和DTMLS 方法的分割結果可能受到損壞，如圖9(f)、圖9(h)、圖9(j)所示.而RECLS 方法對圖9中不同程度灰度不均勻圖像均能保持穩定的分割結果，且分割效率高于LATE 水平集方法以及其它窄帶LATE 方法.

圖9中LATE、DRCTLS、DTMLS 和RECLS 方法對應的迭代次數，運行時間，以及平均每次迭代所需時間如表3所示.DRCTLS 和DTMLS 具有相同的窄帶規模，且DRCTLS 的生成窄帶復雜度高于DTMLS，除圖9(B)的極端情形外，DTMLS 的分割效率總體上高于DRCTLS.

由于DRCTLS 和DTMLS 減少計算區域的收益不足以彌補生成窄帶增加的額外計算開支，反而可能導致總計算效率的下降.RECLS 矩形窄帶生成窄帶復雜度和窄帶面積規模均小于DTMLS 方法，特別是窄帶面積規模的下降，使得窄帶計算效率明顯提升.可見，RECLS 方法的計算效率明顯優于DRCTLS 窄帶法，DTMLS 窄帶法以及未使用窄帶的原始LATE 方法.

3.3 分割精度分析

圖像分割的準確性可以用Jaccard Similarity Coefficient (JSC)[9，10]標準來衡量.

圖9 分割結果對比.第1 行：原圖像；第2 行：LATE 模型分割結果；第3 行：LATE 結合DRCT 窄帶法(DRCTLS)的分割結果；第4 行：LATE 結合DTM 窄帶法(DTMLS)的分割結果；第5 行：我們方法(RECLS)的分割結果.(a)為灰度均勻圖像；(b)-(e)為灰度不均勻圖像；(f)-(j)為灰度嚴重不均勻圖像.

表3 圖9中對應的迭代次數以及運行時間

其中，Ot為標準分割區域，Om為實際分割區域，算符A(·)表示求對應區域的面積.JSC 的取值范圍在0 到1，JSC 的值越大，分割結果越準確.

圖10(a)為二值圖像，圖10(b)-圖10(e)灰度不均勻程度依次遞增.以二值圖10(a)的黑色區域作為標準分割區域Ot，CV 模型，SLGS 模型，RSF 模型，LIC 模型，LSACM 模型，LATE 模型以及本文的RECLS 方法對圖10(a)-圖10(e)的分割結果對應的JSC 值如表4所示.實驗中，圖片尺寸均為100×100.圖10(a)-圖10(c) RECLS 方法的窄帶半徑為5，圖10(d)-圖10(e)RECLS 方法的窄帶半徑為10.

圖10 灰度嚴重不均勻圖像的合成.(a)灰度均勻的二值圖像；(b)(c)灰度不均勻程度一般的合成圖像；(d)(e)灰度嚴重不均勻的合成圖像.

可見，LATE 方法對灰度嚴重不均勻圖像具有較高的分割精度.RECLS 方法對應的JSC 值與LATE 幾乎一致.結合上一節的結論，本文提出的矩形窄帶方法能在不影響LATE 模型分割精度的條件下，提高對灰度嚴重不均勻圖像的分割效率.

4 結論與展望

本文提出一種新的矩形窄帶方法.通過活動約束進一步縮小了窄帶的范圍.利用矩形窄帶代替不規則窄帶，使其更容易與水平集方法相結合，減少了更新水平集的計算量.實驗表明，本文的方法即使在灰度嚴重不均勻情形下也能夠保持穩定的分割結果.

本文的方法在窄帶演化的過程中，可能存在多個矩形窄帶，而這些窄帶沒有實現并行運算.如何讓多個矩形窄帶區域實現并行運算，進一步提高計算的效率，是我們下一步研究的內容.

表4 7 種方法對圖10(a)-圖10(e)的分割結果對應的JSC 值