Shape Analysis课程笔记07-Discrete Surface Curvature
这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要介绍离散网格上计算各种曲率的相关方法。
Applications of Curvature
在上一节中我们介绍了连续曲面上的曲率。具体地，曲面上最重要的两个曲率分别是Gauss曲率和平均曲率，它们与主曲率以及曲面的第二基本形式有着密切的关系。
曲面上的曲率在工程中有非常多的应用。比如说我们可以把Gauss曲率和平均曲率作为曲面形状的描述子、利用曲率来重建曲面、对曲面进行滤波、指导渲染过程、重新网格化(remeshing)等。
Approximating Curvature
Challenge on Meshes
在离散网格上计算曲率的难点在于曲率的定义依赖于曲面上的导数，而在离散网格上每个面片都是一个平面。因此对于离散曲面的曲率我们只能使用一些近似方法来计算。
Local Tensor Methods
早期计算离散曲面上曲率的方法是直接估计网格上的第二基本形式，然后再导出Gauss曲率和平均曲率。Taubin指出在光滑曲面上可以利用一个矩阵来\(M_\mathbf{p}\)来估计曲率：
进一步可以证明\(M_\mathbf{p}\)的特征向量分别为法向\(\mathbf{n}\)以及两个主方向\(\mathbf{t}_1\)和\(\mathbf{t}_2\)，而其特征值与曲面的主曲率密切相关：
对于离散的情况，我们只需要把积分改成求和并替换掉相应的被积项即可：
这里我们简单对离散的过程进行一下推导。对于曲面上的曲线\(\gamma(s)\)，首先根据Taylor展开有：
\[\gamma(s) = \mathbf{p} + \mathbf{T}(0) s + \frac{1}{2} \kappa_{\gamma(0)} \mathbf{N}(0) s^2 + O(s^3)\]上式表明：
\[\Vert \gamma(s) - \mathbf{p} \Vert_2^2 = s^2 + O(s^3)\] \[2 \mathbf{n}_\mathbf{p}^T (\gamma(s) - \mathbf{p}) = \kappa_{\theta} s^2 + O(s^3)\]因此可以得到方向曲率的计算公式：
\[\begin{aligned} \kappa_{\theta} &= \frac{2 \mathbf{n}_\mathbf{p}^T (\gamma(s) - \mathbf{p})}{s^2} + O(s) \\ &= \frac{2 \mathbf{n}_\mathbf{p}^T (\gamma(s) - \mathbf{p})}{\Vert \gamma(s) - \mathbf{p} \Vert_2^2 + O(s^3)} + O(s) \\ &= \frac{2 \mathbf{n}_\mathbf{p}^T (\gamma(s) - \mathbf{p})}{\Vert \gamma(s) - \mathbf{p} \Vert_2^2} + O(s) \end{aligned}\]
Taubin方法的一个缺陷在于它对曲率的估计容易受到局部噪声的影响。由于Taubin方法是基于Taylor展开来推导的，当网格不够稠密时它对曲率的估计往往不够准确。
在一些更现代的方法中，我们会根据实际任务的需求来设计曲率的估计方法。
Direct Approximation of \(\mathbf{II}\)
假设已知切平面上的基向量\(\mathbf{u}\)和\(\mathbf{v}\)，此时曲面的第二基本形式\(\mathbf{II}\)可以表示为一个矩阵。对于切平面上的向量\(\mathbf{w} = c^1 \mathbf{u} + c^2 \mathbf{v}\)，可以利用矩阵形式的\(\mathbf{II}\)来计算\(\mathbf{w}\)方向的shape operator。
因此在离散网格上我们可以构造一个线性系统来求解每个三角形上的\(\mathbf{II}\)：
最后利用切平面的相对旋转和加权平均的方法就可以得到每个顶点上的\(\mathbf{II}\)：
Structure-Preserving Methods
除了上面介绍的方法我们还可以利用曲率的性质来直接推导离散网格上的曲率计算公式。
Structure-Preserving Gaussian Curvature
Gauss-Bonnet定理(Gauss-Bonnet theorem)是Gauss曲率最重要的性质之一，它指出Gauss曲率的积分是一个拓扑不变量。
利用Gauss-Bonnet定理我们可以推导出离散Gauss曲率的表达式。首先定义顶点上的Voronoi cell：
在Voronoi cell上对Gauss曲率进行积分有：
\[\int_V K dA = 2 \pi \chi(V) - \oint_{\partial V} \kappa_g ds\]可以证明在Voronoi cell上欧拉示性数满足\(\chi(V) = 1\)，而测地曲率的积分在每个三角形上等于转角\(\varepsilon_i\)。因此有：
\[\int_V K dA = 2 \pi - \sum_i \varepsilon_i = \int_V K dA = 2 \pi - \sum_i \theta_i\]即Gauss曲率在顶点Voronoi cell上的积分等于\(2\pi\)与顶点内角和的差，称为angle defect：
类似地，我们也可以推导出离散Gauss-Bonnet定理(discrete Gauss-Bonnet theorem)：
Structure-Preserving Mean Curvature
对于曲线来说，\(\kappa \mathbf{N}\)是缩短曲线长度最快的方式。
类似地，mean curvature normal则可以理解为曲面面积的梯度。
对于离散网格，我们可以认为曲面的面积是一个网格到实数的映射：
这样就可以利用mean curvature normal的性质来推导离散曲面上的平均曲率。首先考虑空间中的一个三角形，我们定义一组正交基分别为\(\mathbf{e}\)，\(\mathbf{e}_\perp\)以及\(\mathbf{n}\)。此时\(\mathbf{p}\)点的坐标以及三角形的面积\(A(\mathbf{p})\)可以表示为：
\[\mathbf{p} = p_e \mathbf{e} + p_\perp \mathbf{e}_\perp + p_n \mathbf{n}\] \[A(\mathbf{p}) = \frac{1}{2} b \sqrt{p_n^2 + p_\perp^2}\]对\(A(\mathbf{p})\)进行求导有：
\[\frac{\partial A}{\partial p_e} = 0\] \[\frac{\partial A}{\partial p_n} = \frac{1}{2} b \cdot \frac{p_n}{\sqrt{p_n^2 + p_\perp^2}}\] \[\frac{\partial A}{\partial p_\perp} = \frac{1}{2} b \cdot \frac{p_\perp}{\sqrt{p_n^2 + p_\perp^2}}\]如果我们进一步假设\(p_n=0\)，则\(A(\mathbf{p})\)仅与\(p_\perp\)有关。此时可以把面积\(A(\mathbf{p})\)的梯度表示为：
\[\nabla A(\mathbf{p}) = \frac{1}{2} b \ \mathbf{e}_\perp\]
根据平面几何，我们可以把三角形底边的高向量\(\mathbf{h}\)使用三个顶点位置向量来表示：
\[\frac{b}{h} = \cot \alpha + \cot \beta\] \[\mathbf{h} = \mathbf{p} - \mathbf{p}_0 = \mathbf{p} - \frac{\mathbf{r} \cot \alpha + \mathbf{q} \cot \beta}{\cot \alpha + \cot \beta}\]
把上面的公式结合起来可以得到使用三角形内几何量表示的面积梯度：
\[\begin{aligned} \nabla A(\mathbf{p}) &= \frac{1}{2} b \ \mathbf{e}_\perp = \frac{1}{2} \frac{b}{h} \mathbf{h} \\ &= \frac{1}{2} \bigg( (\mathbf{p} - \mathbf{r}) \cot \alpha + (\mathbf{p} - \mathbf{q}) \cot \beta \bigg) \end{aligned}\]
在三角网格上我们把顶点相邻的所有三角形的面积梯度加起来就得到了曲面的面积梯度：
\[\nabla A(\mathbf{p}) = \frac{1}{2} \sum_j (\cot \alpha_j + \cot \beta_j) (\mathbf{p} - \mathbf{q}_j)\]
因此上面推导的面积梯度即为mean curvature normal在顶点附近的积分。
Alternative Structures
除了上面介绍的方法外我们还可以利用曲率的其它性质来推导离散曲率的表达式。
