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

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

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

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

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

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

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

水平集方法的計(jì)算速度是決定它使用范圍的一個(gè)重要因素，因?yàn)樗鼘⒂绊懙秸麄€(gè)圖像處理過程的進(jìn)度，而初始化步驟是水平集方法中最為耗時(shí)的一個(gè)過程。初始化過程包括對(duì)圖像中所有點(diǎn)的初始距離的計(jì)算和每個(gè)點(diǎn)在初始曲線上的最相近點(diǎn)的確定。即使在使用窄帶方案時(shí)，也要經(jīng)過幾步迭代后才作重新初始化計(jì)算。如果圖像數(shù)據(jù)量很大，整個(gè)計(jì)算時(shí)間會(huì)相當(dāng)可觀，這樣水平集方法就不能用于實(shí)時(shí)計(jì)算了。

2.4.1 快速行進(jìn)算法

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

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

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

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

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

②活躍點(diǎn)（Active），當(dāng)前波前即將經(jīng)過的點(diǎn)，到達(dá)時(shí)間未知；

③遠(yuǎn)點(diǎn)（Faraway），遠(yuǎn)離波前的點(diǎn)，到達(dá)時(shí)間未知。

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

2.4.2 Hermes算法

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

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

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

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

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

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

2.4.3 快速水平集演化算法

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

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

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

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

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

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

STEP6.從窄帶點(diǎn)集中取出速度最快的窄帶點(diǎn)w作為演化種子點(diǎn)，并標(biāo)記w為活動(dòng)點(diǎn)。

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

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

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