推導Discrete Laplace-Beltrami Operator

梯度 - nabla運算元

梯度之前已經寫過，之前給的梯度的求法都是跟可導函式有關，這裡來討論一下三角形網格上的梯度如何定義。

在三角形的三個頂點處有

${\mathbf{x}_i,\mathbf{x}_j,\mathbf{x}_k}$

，我們假設有分段函式f，

$f(\mathbf{x}_i) = f_i,f(\mathbf{x}_j) = f_j,f(\mathbf{x}_k) = f_k$

，

是三⻆角形的⼀次Lagrange插值基函式：

的特點都是在自身i， j， k 之處都等於1，在另外兩個頂點處都為0。

結合重心座標，對於三角形平面上的一個點

$\mathbf{u}$

：

$f(\mathbf{u}) = f_i B_i(\mathbf{u}) + f_j B_j(\mathbf{u}) + f_k B_k(\mathbf{u})\\$

同時：

$B_i(\mathbf{u}) + B_j(\mathbf{u}) + B_k(\mathbf{u}) = 1$

，那麼我們對這個基函式求梯度的話：

$\nabla B_i(\mathbf{u}) + \nabla B_j(\mathbf{u}) + \nabla B_k(\mathbf{u}) = 0$

，所以我們對上面的式子求梯度並代入得：

$\nabla f(\mathbf{u}) = f_i (-\nabla B_j(\mathbf{u}) - \nabla B_k(\mathbf{u})) + f_j \nabla B_j(\mathbf{u}) + f_k \nabla B_k(\mathbf{u})\\$

$\nabla f(\mathbf{u}) = (f_j - f_i) \nabla B_j(\mathbf{u}) + (f_k - f_i) \nabla B_k(\mathbf{u})\\$

同時我們來看這個

，比如對於

，這個變化最大的方向當然就是跟

$(\mathbf{x}_k - \mathbf{x}_j)$

垂直的方向，再看一眼圖：

所以梯度是：

$\nabla B_i(\mathbf{u}) = \frac{(\mathbf{x}_k - \mathbf{x}_j)^{\perp}}{2 A_T}\\$

其中

$(\mathbf{x}_k - \mathbf{x}_j)^{\perp}$

的方向是表示向量在三角形平面逆時針旋轉90°，

表示三角形面積。想法很簡單，也就是

$(\mathbf{x}_k - \mathbf{x}_j) \cdot h/2 = A_T, B_i(\mathbf{x}_i) = 1$

。那麼函式值的變化是（1-0），而梯度方向的變化是（h-0），再加上方向的考量，所以

$\nabla B_i(\mathbf{u}) = \frac{1 -0}{h-0} \mathbf{n}_h = \frac{(\mathbf{x}_k - \mathbf{x}_j)^{\perp}}{2 A_T}$

。

所以：

$\nabla f(\mathbf{u}) = (f_j - f_i) \frac{(\mathbf{x}_i - \mathbf{x}_k)^{\perp}}{2 A_T} + (f_k - f_i) \frac{(\mathbf{x}_j - \mathbf{x}_i)^{\perp}}{2 A_T}\\$

散度 - Laplace運算元

有了上面的梯度計算，我們可以開始推導散度-Laplace運算元在網格上的公式。首先回憶我們的散度的物理意義其實是‘空間中的這個點是否產生液體，其實也就是來看在此點處向內的箭頭比較多還是向外的箭頭比較多’。散度的計算我們用的是 Laplace 運算元，也就是：

$\Delta = \frac{\partial}{\partial x^2} + \frac{\partial}{\partial y^2} + \frac{\partial}{\partial z^2}$

，對於網格上的一點（當然我們這裡的網格都是已經經過處理了，性質非常良好的網格，不會是那種剛剛scan好的，比如有洞啊，奇奇怪怪的網格），使用高斯公式：

$\int_{A_i} div \mathbf{F(u)} dA = \int_{\partial A_i} \mathbf{F(u)} \cdot \mathbf{n(u)} ds\\$

$A_i， \partial A_i$

入上圖所示，這裡的

區域是 mixed Voronoi cell，也就是周邊三角形的外心，當然三角形的外心可能在三角形之外，我們把它clip到三角形的邊上的中點。之所以選擇 mixed Voronoi cell 是因為它產生的error最小（leads to tight error bounds for the discrete operators as shown in ［Meyer et al。 03］）。

同時還值得注意的一點就是，雖然我們這裡貌似把這些三角形都畫在一個平面上，但是我們要知道，這些三角形與三角形之間，當然可能不在同一個平面上。

繼續使用高斯公式：

$\int_{A_i} \Delta f(\mathbf{u}) dA = \int_{A_i} div \nabla f(\mathbf{u}) dA = \int_{\partial A_i} \nabla f(\mathbf{u}) \cdot \mathbf{n(u)} ds\\$

$\nabla f$

在每個三角形上是常量，所以我們對單個三角形來計算：

$\int_{\partial A_i \cap T} \nabla f(\mathbf{u}) \cdot \mathbf{n(u)} ds = \nabla f \cdot (\mathbf{a} - \mathbf{b})^{\perp} = \frac{1}{2} \nabla f \cdot (\mathbf{x}_j - \mathbf{x}_k)^{\perp}\\$

這裡是因為我們選取的是三角形的外心，所以當然

$\mathbf{a}, \mathbf{b} , (\mathbf{x}_j - \mathbf{x}_k)$

有以上關係。

我們繼續代入上面求得的關於梯度的表示式：

$\int_{\partial A_i \cap T} \nabla f(\mathbf{u}) \cdot \mathbf{n(u)} ds = (f_j - f_i) \frac{(\mathbf{x}_i - \mathbf{x}_k)^{\perp} \cdot (\mathbf{x}_j - \mathbf{x}_k)^{\perp} }{4 A_T} + (f_k - f_i) \frac{(\mathbf{x}_j - \mathbf{x}_i)^{\perp} \cdot (\mathbf{x}_j - \mathbf{x}_k)^{\perp} }{4 A_T}\\$

用

$\gamma_j, \gamma_k$

來表示

處的角度，我們有：

$A_T = \frac{1}{2} \sin \gamma_j |\mathbf{x}_j - \mathbf{x}_i| |\mathbf{x}_j-\mathbf{x}_k| = \frac{1}{2} \sin \gamma_k |\mathbf{x}_i - \mathbf{x}_k| |\mathbf{x}_j-\mathbf{x}_k|\\$

又根據點乘關係：

$\cos \gamma_j = \frac{(\mathbf{x}_j - \mathbf{x}_i) \cdot (\mathbf{x}_j - \mathbf{x}_k)}{|\mathbf{x}_j - \mathbf{x}_i| |\mathbf{x}_j - \mathbf{x}_k|}\\$

$\cos \gamma_k = \frac{(\mathbf{x}_i - \mathbf{x}_k) \cdot (\mathbf{x}_j - \mathbf{x}_k)}{|\mathbf{x}_i - \mathbf{x}_k| |\mathbf{x}_j - \mathbf{x}_k|}\\$

所以上面的式子可以化簡為：

$\begin{align} \int_{\partial A_i \cap T} \nabla f(\mathbf{u}) \cdot \mathbf{n(u)} ds &= (f_j - f_i) \frac{\cos \gamma_k |\mathbf{x}_i - \mathbf{x}_k| |\mathbf{x}_j - \mathbf{x}_k| }{2 \sin \gamma_k |\mathbf{x}_i - \mathbf{x}_k| |\mathbf{x}_j-\mathbf{x}_k| } + (f_k - f_i) \frac{ \cos \gamma_j|\mathbf{x}_j - \mathbf{x}_i| |\mathbf{x}_j - \mathbf{x}_k|}{2 \sin \gamma_j |\mathbf{x}_j - \mathbf{x}_i| |\mathbf{x}_j-\mathbf{x}_k|} \\ &= \frac{1}{2} (\cot \gamma_k (f_j - f_i) + \cot \gamma_j (f_k - f_i)) \end{align} \\$