2.4 Allen-Cahn去噪模型的水平集求解

水平集方法的主要思想是把二維平面的閉合曲線視為三維空間連續(xù)函數(shù)曲面φ（水平集函數(shù)）的一個水平層，通常是{φ=0}的一個水平層。這樣，二維曲線的演化就轉(zhuǎn)換為三維水平集函數(shù)曲面的演化。要使得水平集函數(shù)曲面的演化結(jié)果和閉合曲線的演化相關(guān)，水平集函數(shù)曲面的演化要遵循如下Hamilton-Jacobi類型的偏微分方程：

式中，F(xiàn)是和曲線在其法線方向的速度相關(guān)的某種數(shù)量，通常情況下F只在曲線的位置上有意義，也就是定義在零水平層上，因此需要將零水平層的F擴展到整個函數(shù)曲面定義域，是擴展速度場Fext。一般來說，有兩種速度擴展的方法。一種是全局擴展算法，另一種是窄帶算法，速度比全局算法快。這兩種算法的過程基本上是一樣的，只是要處理的數(shù)據(jù)量不一樣。窄帶算法需要增加與限制數(shù)據(jù)量的相關(guān)操作。窄帶法最初由Chop提出，Adalsteinsson給出了詳細的實現(xiàn)方法。基本思想是只演化位于零水平集周圍很窄的一個帶狀區(qū)域水平集函數(shù)的值。基本的過程描述如下：

①初始化平面上的所有點，包括計算初始距離和執(zhí)行速度擴展；

②用速度函數(shù)作迭代計算；

③跟蹤曲線經(jīng)過一步迭代后的新位置；

④檢查終止條件是否滿足，如果不滿足則進行新一次的迭代計算。

水平集方法的計算速度是決定它使用范圍的一個重要因素，因為它將影響到整個圖像處理過程的進度，而初始化步驟是水平集方法中最為耗時的一個過程。初始化過程包括對圖像中所有點的初始距離的計算和每個點在初始曲線上的最相近點的確定。即使在使用窄帶方案時，也要經(jīng)過幾步迭代后才作重新初始化計算。如果圖像數(shù)據(jù)量很大，整個計算時間會相當可觀，這樣水平集方法就不能用于實時計算了。

2.4.1 快速行進算法

在窄帶算法的計算過程中，第一步就是初始化步驟。一種最簡單的初始化方法就是通過計算當前點與曲線上所有點的距離，尋找距離的最小值，以得到與當前點距離最近的曲線上的點，同時記錄下最近點的坐標。這種方法稱為直接法。直接法的計算速度為O（N×n），這里N是平面上點的數(shù)量，n是曲線上點的數(shù)量。Adalsteinsson提出了用快速行進法為水平集作初始化計算，這是一種用來求解Eikonal方程的數(shù)值算法。算法中結(jié)合了Hamilton-Jacobi方程的黏性解法中的上推算法、水平集方法中的窄帶算法和一個最小數(shù)據(jù)堆的數(shù)據(jù)結(jié)構(gòu)。

設(shè)C是一個初始曲面，F(xiàn)是沿曲面法線方向的速度。考慮當曲線經(jīng)過點（x，y，z）的特殊情況，這時函數(shù)T（x，y，z）滿足方程|?F|T=1。根據(jù)Osher-Sethian求解水平集的“熵守桓”差分方法，可得方程|?F|T=1的數(shù)值解表達式如下：

式中，D-和D+分別是前向差分算子和后向差分算子。由于上式從本質(zhì)上說是每個點的二次方程，因此可以用迭代的方法求解方程，然后選擇最大可能的解作為方程最后的解，這和水平集方程的黏性解法是一致的。

由式（2.12）可知，邊界傳播方向是從T比較小的點流向T比較大的點，根據(jù)這樣的特點，Sethian提出了快速行進算法來迅速傳播邊界值T。基本思想是在傳播邊界外圍構(gòu)造一個激活（Active）窄帶，窄帶內(nèi)的點的到達時間未定，利用迎風機制（Upwind Scheme）將當前邊界向外傳播，就像水波擴散一樣，凡是擴散到的點，就凍結(jié)其波前到達時間，然后再根據(jù)當前的波前構(gòu)造新的激活帶，如此循環(huán)，就可以得到整個平面上每點的到達時間。因此，建立快速行進算法的關(guān)鍵是：用上推的方式掃描當前曲線周圍一個窄帶范圍的點，得到曲線演化最先到達的點，從而把曲線不斷地向前推進。根據(jù)計算網(wǎng)格每點距離閉合曲線的遠近，將網(wǎng)格點分為三類：

①激活點（Alive），已經(jīng)被波前傳播過的點，到達時間已知；

②活躍點（Active），當前波前即將經(jīng)過的點，到達時間未知；

③遠點（Faraway），遠離波前的點，到達時間未知。

快速算法的主要技術(shù)是一個最小堆的數(shù)據(jù)結(jié)構(gòu)。這個結(jié)構(gòu)實際上是一個完全二叉樹，二叉樹中每個父節(jié)點的值要不大于其子節(jié)點的值。在實際應(yīng)用中，一般會把二叉樹存儲在一個數(shù)組中，而子節(jié)點的位置正好是父節(jié)點的兩倍。因此，存儲最小值的根節(jié)點是數(shù)組的第一個元素。尋找父節(jié)點或子節(jié)點的位置只要花費O（1）的時間。在存儲時間T的同時，需要存儲點在網(wǎng)格中的位置坐標。快速行進算法通過不斷尋找窄帶中的最小值（FindSmallest）而達到求解的目的。FindSmallest操作包括去掉根結(jié)點和一次對二叉樹的下推（DownHeap）掃描，DownHeap操作是為了保證剩余的點滿足二叉樹的性質(zhì)。算法步驟為：標記當前點的鄰點中當前狀態(tài)是FarAway的狀態(tài)為Trial，這是通過插入（Insert）操作完成的；同時還要重新計算當前點的鄰點的到達時間。Insert增加節(jié)點的數(shù)量，并且用上推（UpHeap）操作把插入點放在合適的位置以保證二叉樹的性質(zhì)。快速行進算法中的很重要一點在于需要在計算過程中始終保持二叉樹組的性質(zhì)。這可以使我們在查找窄帶中最小值點的過程中值只需要O（1）的時間就可以完成。這樣，遍歷整個圖像點的過程需要花費O（Nlog（N））的時間。

2.4.2 Hermes算法

FMM算法是一種有效降低水平集計算開銷的方法，它通過計算Eikonal方程得到網(wǎng)格點的行進速度，通過這種方式可以使計算復(fù)雜度從O（N2）降至O（Nk），其中N為水平（或垂直）方向的網(wǎng)格點數(shù)目，k為窄帶的寬度，但該算法要求演化種子點的行進方向保持不變。Paragios提出的Hermes算法對經(jīng)典FMM算法進行改進，在保留了FMM算法快速優(yōu)點的同時，解除了行進速度方向保持不變的限制。

本章將速度函數(shù)應(yīng)用到Hermes算法中，設(shè)計出式（2.10）的水平集方程數(shù)值算法。

式（2.10）在Hermes算法中的離散差分為

式中，h和τ分別為空間步長和時間步長，V為演化曲線的速度函數(shù)。對流項?g·?u的迎風型有限差分逼近和Δu項的拉普拉斯差分逼近為

式中，，D+和D-分別表示前向差分和后向差分算子。

行進速度函數(shù)V的有限差分表示為

2.4.3 快速水平集演化算法

根據(jù)Hermes算法和演化曲線速度函數(shù)V的描述，設(shè)計APNACM算法如下。

STEP1.輸入含噪圖像，設(shè)定移動界面寬度參數(shù)ξ和積分閾值參數(shù)γ。

STEP2.含噪輸入圖像中的所有網(wǎng)格點標記為遠離點，并將速度值v置為∞。

STEP3.計算初始圖像網(wǎng)格點（ih，jh）的u值，確定初始活動點集。如果網(wǎng)格點（ih，jh）的u≤0.5，且與其相鄰的4個網(wǎng)格點中的任何一個u＞0.5，則認為該點是初始圖像曲線所在網(wǎng)格點，并將該點標記為活動點。

STEP4.提取STEP3中生成的一個活動點v，根據(jù)式（2.16）依次對v相鄰的4個網(wǎng)格點wi（i=1，2，3，4）計算wi的速度值vi：如果wi是窄帶點，則更新窄帶點wi的速度值；如果wi不是窄帶點，則標記該點為窄帶點，將其插入窄帶點集。

STEP5.如果STEP3中生成的所有活動點都處理完畢，則執(zhí)行STEP6；否則返回STEP3。

STEP6.從窄帶點集中取出速度最快的窄帶點w作為演化種子點，并標記w為活動點。

STEP7.如果水平集曲線積分值小于閾值參數(shù)γ，則轉(zhuǎn)至STEP9，否則依次對w點相鄰的4個網(wǎng)格點wi（i=1，2，3，4）進行標記和行進速度的更新：如果①wi是活動點，則返回STEP7；②wi是窄帶點，則根據(jù)式（2.16）計算wi的新速度值vi，同時更新窄帶點集中wi的相應(yīng)速度值vi；③wi不是窄帶點，則根據(jù)式（2.16）計算wi的新速度值vi，同時標記wi為窄帶點，將其插入窄帶點集。

STEP8.如果窄帶點集不空，則返回STEP6，繼續(xù)下一個窄帶點演化；否則執(zhí)行STEP9。

STEP9.輸出去噪圖像，結(jié)束。