はじめに

三角形要素で分割された曲面上での曲率の計算をしたく、下記の論文の手法を解説するのが本記事の趣旨です。
Meyer M., Desbrun M., Schröder P., Barr A.H. (2003) Discrete Differential-Geometry Operators for Triangulated 2-Manifolds. In: Hege HC., Polthier K. (eds) Visualization and Mathematics III. Mathematics and Visualization. Springer, Berlin, Heidelberg
https://people.eecs.berkeley.edu/~jrs/meshpapers/MeyerDesbrunSchroderBarr.pdf

図が無いとわかりにくいので、そのうち追加します。みんなこの手の絵ってどうやって書いてるんだろう…。

この方法でのガウス分布の実装については、下の記事の方が手っ取り早いと思います。
http://mikuhatsune.hatenadiary.com/entries/2016/04/30

前提条件

三次元空間上にある曲面$S$を考え，曲面を三角形一次要素（3節点）で有限要素分割する。格子の各頂点（節点）の座標を$\mathbf{x}\in\mathbb{R}^3$とし、回転自由度などは与えない。節点に回転自由度を与えれば、各格子ごとに曲率を陽に定義できるが(discrete Kirchhoff triangle(DKT)要素など)、ここではこの手の方法は考慮しない。加えて、各頂点について、周囲の頂点の空間分布に基づき二次曲面を最小二乗法などで近似する方法も用いない。あくまで格子分割された情報を使う。

定義

一般的な微分幾何の入門書の範囲で説明される定義。自身が無い人は下記などに目を通しておくと直感的にわかりやすいと思います。
https://www.kyoritsu-pub.co.jp/bookdetail/9784320017887

曲面内の任意の位置で法線ベクトル$\bf{n}$・法線ベクトルと直交する平面（接平面）を定義する。接平面上の任意の単位方向ベクトルを$\mathbf{e}_{\theta}$とする。
$\mathbf{n}$と$\mathbf{e}_{\theta}$で定義できる平面と$S$が重なってできる曲線の曲率を法曲率$\kappa^{N}(\theta)$と定義する。
主曲率$\kappa_1, \kappa_2$と、対応する単位方向ベクトル$\mathbf{e}_1, \mathbf{e}_2$を定義する。

法曲率$\kappa^N$と主曲率との関係は下記の通り（法曲率のオイラーの公式）。　

\kappa_N(\theta)=\kappa_1\cos{\theta}^2+\kappa_2\sin{\theta}^2.

証明は下記参照。
http://sshmathgeom.private.coocan.jp/diffgeom/eulerdupin.pdf

この関係を使うと，平均曲率$\kappa_H$は下記のように定義される。

\kappa_H=\frac{1}{2\pi}\int_0^{2\pi}{\kappa^N(\theta)d\theta}=\frac{\kappa_1+\kappa_2}{2}.

最後に、Gauss曲率$\kappa_G$の定義は下記の通り。

\kappa_G=\kappa_1\kappa_2.

Laplace-Beltrami作用素

Laplace-Beltrami作用素はLaplace作用素の一般座標系への拡張であり、本気で考えるとリーマン計量の世界に入らざるを得ない。よって、この節も定義だけを並べておく。

Laplace-Beltrami作用素を$\Delta_S$とすると、曲面上の座標$\mathbf{x}_P$において、

\Delta_s\mathbf{x}_P=-2\kappa_H\mathbf{n}=\mathbf{K}(\mathbf{x}_P)

が成立する。後の記述の簡略化のため、$\mathbf{K}$を導入した。

座標$\mathbf{x}_P$において，点の周囲の微小面積$\mathcal{A}$を与えると，下記の解釈も可能となる。

2\kappa_H\mathbf{n}=\lim_{diam(\mathcal{A}))\rightarrow 0}\frac{\nabla\mathcal{A}}{\mathcal{A}}

ここで$diam(\mathcal{A})$は直径を表す．

#　空間離散化
各節点について、節点を共有する三角形要素から適切に表面積を選択する。対象とする節点を端点に持つ辺の中点と、三角形要素の重心もしくは外心を直線で結ぶことでできる閉曲面（1-ring）を用いることとし、その内部を対象とする表面積と定義する。重心か外心かはどちらを選んでもよいが、三角形の角度がいずれも鋭角であれば、外心を用いた方が精度が良いことが論文中で示されている。ここでは記述の簡略化のために下記以降は三角形要素一つの面積を$\mathcal{A}_T$と表す。

上で定義した面積を利用して、$\mathbf{K}$の面積平均を離散的に計算する。この演算は

\begin{align}
\int_{\mathcal{A}}{\mathbf{K}dA}&=-\int_{\mathcal{A}}{\Delta_S\mathbf{x}dA} \\
&=-\int_{\mathcal{\partial}\mathcal{A}}{\mathbf{n}\cdot\nabla_S\mathbf{x}d\Gamma}
\end{align}

として、境界積分の形に書き直せる（Gauss-Greenの定理）。この演算を対象となる三角形ごとに離散的に行い、後で足し合わせる。法線ベクトルの境界積分について、対象とする三角形Tの各節点の座標を$\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3$とし、今回対象とする節点を$\mathbf{x}_1$とすれば、計算に用いる1-ring内の辺は節点2-3の辺のみ。そこで、

\begin{align}
\int_{\partial\mathcal{A}_M\cap T}{\mathbf{n}\cdot\nabla_S\mathbf{x}d\Gamma}&=\frac{1}{2}[\mathbf{n}\times(\mathbf{x}_2-\mathbf{x}_3)]\cdot\nabla_S\mathbf{x} \\
&=\frac{1}{2}[\mathbf{x}_2-\mathbf{x}_3]^{\perp}\cdot\nabla_S\mathbf{x}
\end{align}

と書ける。ここで$\perp$は反時計回りに90°回転させる演算を表す。これは幾何学的状況を考えると明らか。

今回は三角形一次要素を対象としているため、$\nabla_S\mathbf{x}$は一定値であり，積分区間から外せる。$\nabla_S\mathbf{x}$の計算には有限要素法における三角形一次要素の形状関数$N_i$を導入する。三次元空間において、三角形内部に局所座標系を定義し、勾配を整理すれば

\begin{align}
\nabla_S N_i&=\frac{1}{2\mathcal{A}_T}[(\mathbf{n}\times(\mathbf{x}_k-\mathbf{x}_j)] \\
&=\frac{1}{2\mathcal{A}_T}[\mathbf{x}_k-\mathbf{x}_j]^{\perp}
\end{align}

となる($i\rightarrow j\rightarrow k$の順)。通常の有限要素法の計算ではここまで書き下すことは少ないので、少し計算が煩雑かもしれない。

形状関数の定義から

\nabla_S N_1+\nabla_S N_2+\nabla_S N_3=\mathbf{0}

が自明なので、

\begin{align}
\nabla_S\mathbf{x}&=\nabla_S N_2(\mathbf{x}_2-\mathbf{x}_1)+\nabla_S N_3(\mathbf{x}_3-\mathbf{x}_1) \\
&=\frac{1}{2\mathcal{A}_T}[\mathbf{x}_1-\mathbf{x}_3]^{\perp}(\mathbf{x}_2-\mathbf{x}_1)^T+\frac{1}{2\mathcal{A}_T}[\mathbf{x}_2-\mathbf{x}_1]^{\perp}(\mathbf{x}_3-\mathbf{x}_1)^T
\end{align}

とする。これを先程の面積積分に代入すれば、

\begin{align}
\int_{\partial\mathcal{A}_M\cap T}{\mathbf{n}\cdot\nabla_S\mathbf{x}d\Gamma}&= \\
\frac{1}{4\mathcal{A}_T}&\left[[\mathbf{x}_2-\mathbf{x}_3]^{\perp}\cdot[\mathbf{x}_1-\mathbf{x}_3]^{\perp}(\mathbf{x}_2-\mathbf{x}_1)+[\mathbf{x}_2-\mathbf{x}_3]^{\perp}\cdot[\mathbf{x}_2-\mathbf{x}_1]^{\perp}(\mathbf{x}_3-\mathbf{x}_1) \right]\\
\end{align}

である。このとき

[\mathbf{x}_2-\mathbf{x}_3]^{\perp}\cdot[\mathbf{x}_1-\mathbf{x}_3]^{\perp}=\|\mathbf{x}_2-\mathbf{x}_3\|\|\mathbf{x}_1-\mathbf{x}_3\|\cos{\theta_{231}},

かつ

\begin{align}
\mathcal{A}_T&=\frac{1}{2}\|(\mathbf{x}_2-\mathbf{x}_3)\times(\mathbf{x}_1-\mathbf{x}_3)\| \\
&=\frac{1}{2}\|\mathbf{x}_2-\mathbf{x}_3\|\|\mathbf{x}_1-\mathbf{x}_3\|\sin{\theta_{231}}
\end{align}

だから、計算を大幅に簡略化でき、整理すれば

\begin{align}
\int_{\partial\mathcal{A}_M\cap T}{\mathbf{n}\cdot\nabla_S\mathbf{x}d\Gamma}= 
\frac{1}{2}\left[\cot{\theta_{231}}(\mathbf{x}_2-\mathbf{x}_1)+\cot{\theta_{123}}(\mathbf{x}_3-\mathbf{x}_1) \right]\\
\end{align}

となる。この式を全ての三角形で足し合わせると、

\begin{align}
\int_{\mathcal{A}_M}{\mathbf{K}dA}&=\frac{1}{2}\sum_{j\in N_1(i)}{\left[\cot{\theta_{+}}+\cot{\theta_{-}}\right](\mathbf{x}_i-\mathbf{x}_j)}
\end{align}

となる。右辺は評価点を含む全ての辺に対する総和の形で表せ、$\theta_{+}$と$\theta_{-}$は辺を挟むそれぞれの三角形が作る角度を表す。

離散平均曲率・法線

定義より、

\mathbf{K}=\frac{1}{2\mathcal{A}_{all}}\sum_{j\in N_1(i)}{\left[\cot{\theta_{+}}+\cot{\theta_{-}}\right](\mathbf{x}_i-\mathbf{x}_j)}

として離散曲率法線を求め、このベクトルの絶対値の半分を平均曲率とする。また、単位ベクトルをとればそれが法線ベクトルに対応する。

上の方法の場合、平均曲率はベクトルの絶対値の形で出力されるため符号を持たない。そもそも法線ベクトルの計算において三角形要素の辺の情報を用いており、先見的に与えられる三角形要素の外向き法線ベクトルの情報を用いていない。上の計算後、法線ベクトルの向きを補正した後、正負を入れ替える必要のある節点の曲率に負の符号を設定する。

計算例

上の式を実装して計算してみました（下図）。定性的には妥当そうです。

定式化を見るに、三角形要素の面積計算などでざっくりした計算を導入しているし、定量性は微妙そうな…。

備考

三角形要素の平均曲率法線
https://u-1roh.hatenadiary.org/entry/20060711/1152627733
http://tkenichi.hatenablog.jp/entry/2015/05/23/180109