如何計算圖像的曲率？

01-26

一般的數字圖像都是一個二維離散函數 $I(x_i,y_j)$ ，圖像處理的任務通常是想得到一個新的圖像 $U(x_i,y_j)$ 。這個 $U(x_i,y_j)$ 滿足特定的性質，比方說去模糊、去霧、超解析度、去噪、分割等等。以去模糊為例，得到的圖像要比原始圖像清晰。這個過程我們通常要建立一個數學模型 ${cal E}(U)$ ，使得該模型最小值對應的 $U(x_i,y_j)$ 即為所求。顯然，模型的建立對於最終的結果至關重要。另一方面，模型也要易於求解才行。所以說，模型和求解方法是一對孿生姐妹。

最為通用的求解方法是梯度下降法，即: $frac{partial U^t}{partial t}=-frac{partial {cal E}}{partial U}$

但梯度下降法也是最為低效率的方法之一。原因是多方面的。首先，它沒有考慮 ${cal E}(U)$ 的可導性。對於不可導的模型可能需要做光滑逼近。其次，在時間步長的選取上它需要滿足數值穩定性要求，即便是選取高階精度的時間離散化公式（梯形公式、龍哥庫塔方法等等）。再次，它沒有考慮模型 ${cal E}(U)$ 的具體形式，所以雖然非常通用，但對於具體問題又非常沒有效率。

另外一種求解方法是從如下必要條件構造迭代公式： $frac{partial {cal E}}{partial U}=0$

通常，這種方法比梯度下降法更有效率，也需要了解模型 ${cal E}(U)$ 的具體形式。但構造這種方法一般比較困難，得到的結果很難看出其對應的物理意義。

前文基於曲率的圖像處理介紹的曲率濾波方法針對曲率正則項設計了一種全新的優化方法，能夠高效地減少圖像的曲率而不需要計算曲率。詳細的介紹見曲率濾波博士論文第六章。針對圖像的曲率優化問題，這一演算法是目前最高效的優化方法（沒有之一）。

那麼，計算圖像的曲率是不是就不必要了呢？其實也不是。最明顯的例子還是曲率濾波。雖然我們能優化曲率而不需要計算曲率，我們還是需要觀察模型 ${cal E}(U)$ 以確保收斂。這就需要計算曲率正則項的能量。如何計算曲率呢，尤其是針對離散的數字圖像？

讓我們以平均曲率為例。我們考慮曲面 $Psi(vec{x})=(vec{x},U(vec{x}))$ 。我們很容易計算它的第一和第二基本形式，從而得到它的平均曲率公式：

$H=frac{(1+U_x^2)U_{yy}-2U_xU_yU_{xy} +(1+U_y^2)U_{xx}}{2(1+U_x^2+U_y^2)^{frac{3}{2}}}$

大部分跟平均曲率相關的論文都通過離散化這個公式來計算平均曲率。另外一種計算方法是通過二次曲面擬合來避免離散化，即

$U(x,y)approx f(x,y)equiv C_5x^2+C_4y^2+C_3xy+C_2x+C_1y+ C_0$

然後通過最小二乘法先確定係數 $C_i$ ，之後將 $C_i$ 代入上式，得到：

$Happrox frac{(1+C_2^2)C_4-C_2C_1C_3 +(1+C_1^2)C_5}{(1+C_1^2+C_2^2)^{frac{3}{2}}}$

這種方法的計算量比較大，而且沒有用到圖像regular sampling的性質，所以在圖像處理中用得不多，而是在三維點雲和三角網格等不規則採樣的處理中比較常見。

讀到這裡，聰明的讀者會問：有沒有不需要離散化而且利用規則採樣的計算公式呢？筆者在構造平均曲率濾波的過程中，不經意地發現了一種全新的計算公式（曲率濾波博士論文的公式6.12）。希望有經驗的讀者斧正：

$Happrox left( egin{array}{ccc} frac{-1}{16} & frac{5}{16} &frac{-1}{16}\ frac{5}{16} & -1 &frac{5}{16}\ frac{-1}{16} & frac{5}{16} &frac{-1}{16}\ end{array} ight)ast U$

這個公式是一個線性卷積操作，跟上面兩個非線性的計算公式完全不同。這個公式是如何構造出來的呢？

在解釋這個公式之前，讓我們回憶一下微分幾何中大名鼎鼎的歐拉理論（Euler Theorem, 1760，維基百科，好吧其實我當初看的是陳省身的《微分幾何初步》）： $kappa_{ heta}=kappa_1cos^2 heta+kappa_2sin^2 heta$ ，

其中 $kappa_{ heta}$ 是方向曲率， $kappa_1$ 和 $kappa_2$ 為主曲率， $heta$ 是到主平面的夾角。我們可以對 $heta$ 求積分，得到

$egin{split}intlimits_{-pi}^{pi}kappa_{ heta}mathrm{d} heta=&kappa_1intlimits_{-pi}^{pi}cos^2 hetamathrm{d} heta+kappa_2intlimits_{-pi}^{pi}sin^2 hetamathrm{d} heta\=&kappa_1pi+kappa_2pi=pi(kappa_1+kappa_2)=2pi Hend{split},.$

從而，我們有：

$H=frac{1}{2pi}intlimits_{-pi}^{pi}kappa_{ heta}mathrm{d} heta$

注意，這個公式是數學上嚴格成立的，所以不是約等於符號，是嚴格等於。這個公式說，平均曲率是所有方向曲率的平均。這也是它為什麼被稱為「平均」曲率的原因（不光指它是兩個主曲率的平均）。筆者在曲率濾波博士論文中計算各個方向上的方向曲率 $kappa_{ heta}$ ，然後相加，再除以 $2pi$ ，最終得到上述線性卷積公式。顯然，方向曲率是被離散逼近的，所以上述線性計算公式是約等於符號。這裡我們利用了局部窗口中方向曲率是有限的、離散的這一特性。我們的計算公式是線性的，充分利用了方向的離散型，從而避免了使用連續公式然後離散化帶來的誤差（計算是直接從離散圖像來的）。當然，曲率濾波博士論文中方向曲率的計算是一個線性近似公式，然後假設近似誤差在多次平均的情況下相互抵消（如果你想到了大數定理，說明你的數學比我好）。感興趣的讀者可以看看曲率濾波博士論文第六章。

那麼線性公式和非線性公式哪個的計算精度更高呢？筆者在曲率濾波博士論文中做了簡單的測試，用兩個公式分別計算一個球圖像的曲率，然後跟真實值比較。得到的結論是線性公式更為準確。當然，這個實驗也只是針對球的圖像，並不是所有圖像都會有相同的結論。不過我們在500幅自然圖像（BSDS500）上的對比結果表明，兩者估計的平均曲率差別不大。但是從計算量和公式的優美程度上來講，筆者還是傾向於線性卷積公式。

上圖左邊為原始圖像，中間為非線性經典公式計算的結果，右邊為線性卷積公式計算的結果。

線性卷積公式跟經典的非線性公式相比，另外一個優勢是它對 $U(x_i,y_j)$ 完全沒有任何要求，而非線性公式則要求圖像至少二次可導地光滑。這也導致非線性公式在大梯度附近不能正確地估計平均曲率。顯然，沒有光滑性要求的線性公式更適合估計平均曲率，尤其是在大梯度區域。

一個是非線性的經典公式，一個是線性的卷積公式，兩者都可以用來估計平均曲率。你會選擇哪一個呢？

最後，再次對比兩個公式：

$H=frac{(1+U_x^2)U_{yy}-2U_xU_yU_{xy} +(1+U_y^2)U_{xx}}{2(1+U_x^2+U_y^2)^{frac{3}{2}}}$ ，
$Happrox left( egin{array}{ccc} frac{-1}{16} & frac{5}{16} &frac{-1}{16}\ frac{5}{16} & -1 &frac{5}{16}\ frac{-1}{16} & frac{5}{16} &frac{-1}{16}\ end{array} ight)ast U$ （曲率濾波博士論文中的公式6.12）

丙申年秋

參考文獻：

@phdthesis{gong:phd, title={Spectrally regularized surfaces}, author={Gong, Yuanhao}, year={2015}, school={ETH Zurich, Nr. 22616}, note={http://dx.doi.org/10.3929/ethz-a-010438292}}@article{gong:Bernstein, Author = {Yuanhao Gong}, Journal = {ICASSP}, Month = {March}, Pages = {1701--1705}, Title = {Bernstein Filter: a new solver for mean curvature regularized models}, Year = {2016}}