多边形网格处理（4）

4. Smoothing（平滑）

• 网格平滑（mesh smoothing）
  • 从抽象的层面看，网格平滑是指设计和计算一个三角形的光滑函数$\mathbf{f}:\mathcal{S}\to {\mathbb{R}}^d$
  • Mesh smoothing 是几何处理的基础工具
  • 光滑函数可以使用，例如顶点位置、纹理坐标、或顶点偏移量来描述
  • 可用于网格参数化、各向异性remeshing、空洞填补、变形等任务
• 网格平滑的两块内容：去噪和光顺
  • 去噪的意思是从光滑函数$\mathbf{f}$中去除高频噪声。大多数情形下，$\mathbf{f}$指的是顶点位置，通常会被高频噪声影响。去噪要求就是去除噪声（高频部分），同时保持整体形状（低频部分），需要拓展频率和低通滤波器在离散三角网格上的概念
  • 光顺，不单是平滑函数$\mathbf{f}$，还要获得尽可能平滑的表面块（as-smooth-as-possible surface patch），或尽可能平滑的形状变换（as-smooth-as-possible shape deformation）。as-smooth-as-possible 的意思就是某种光顺能量的最小化，比如考虑了曲率或高阶微分量的

4.1 频域分析工具在离散三角网格上的拓展

4.1.1 1D Fourier Transform

• 傅里叶变换把一个单变量函数从空域变换到频率域：

$$F(\omega )={\int }_{-\infty }^{\infty }f(x)e^{-2\pi i\omega x}dx$$

$$f(x)={\int }_{-\infty }^{\infty }F(\omega )e^{2\pi i\omega x}d\omega$$

这些方程有一种直观的几何解释：

• 函数$f(x)$可以看做是某个特定向量空间的元素，可以使用内积：

$$<f,g>={\int }_{-\infty }^{\infty }f(x)\overline{g(x)}dx$$

• 复指数函数可以看做是频率为$\omega$的复波，它们构成了频率相关的正交基

$$e_{\omega }(x):=e^{2\pi i\omega x}=\cos (2\pi \omega x)-i\sin (2\pi \omega x)$$

• 在这种情况下，傅里叶变换可以看做是换基，通过正交投影把向量$f$投影到到基向量们$e_{\omega }$上：

$$f(x)=\sum _{\omega =-\infty }^{\infty }<f,e_{\omega }>e_{\omega }$$

尺度系数$<f,e_{\omega }>$就是$F(\omega )$，表示$f$中包含了多少关于基函数$e_{\omega }$的成分，也就是信号$f(x)$中关于频率$\omega$的振幅是多少。因为$\omega$是连续的，所以求和进一步变成积分，就是式（2）了

• 因为$F(\omega )$所在坐标系直接反映了频率，因此简单的提出超过某个阈值的频率，就是理想低通滤波器，这与重建滤波后的函数等价：

$$\tilde{f}(x)={\int }_{-{\omega }_{max}}^{{\omega }_{max}}<f,e_{\omega }>e_{\omega }d\omega$$

4.1.2 流形谐波

现在把1D Fourier框架拓展到二维流形上。傅里叶变换不能直接转换到流形函数上。但可以观察到：

• 正弦和余弦函数，也就是复波$e_{\omega }$，是Laplace算子的特征方程：

$$\Delta (e^{2\pi i\omega x})=\frac{d^2}{d{x}^2}{e}^{2\pi i\omega x}=-(2\pi \omega {)}^2{e}^{2\pi i\omega x}$$

一个方程$e_i$是Laplacian的特征方程，若$\Delta e_i={\lambda }_i{e}_i$

• 由此可见，1D Fourier transform的基函数是Laplacian的特征方程。因此，很自然地选择Laplace-Beltrami算子的特征方程作为泛化的基函数。已知如何离散化Laplace-Beltrami，也就提供了在离散三角网格上泛化Fourier transform的可能
• 这个思路同样可用于其他双变量基函数：对2D矩形，Laplace-Beltrami算子的特征方程对应了离散余弦变换的基函数；对球面，对应了球谐波
• 因此，Laplace-Beltrami的特征方程扩展到了任意二维流形上，称做流形谐波（manifold harmonics）

对离散化的三角网格，将连续函数$f(\mathbf{x})$表示为网格顶点的采样值组成的向量$(f(v_1),\cdots ,f(v_n){)}^T$，有
$$\left(\begin{array}{c}\Delta f(v_1)\\ ...\\ \Delta f(v_n)\end{array}\right)=\mathbf{L}\left(\begin{array}{c}f(v_1)\\ ...\\ f(v_n)\end{array}\right)$$

• Laplace-Beltrami matrix，其中每一行的权重是顶点$v_i$的离散Laplacian：

$$\Delta f(v_i)=\sum _{v_j\in {\mathcal{N}}_1(v_i)}{w}_{ij}(f(v_j)-f(v_i))$$

特征方程$e_{\omega }(x)$变成了Laplace矩阵的特征向量${\mathbf{e}}_1,\cdots ,{\mathbf{e}}_n$：一个n为特征向量${\mathbf{e}}_i$可以看做连续特征方程${\mathbf{e}}_i(\mathbf{x})$的离散采样$({\mathbf{e}}_i(v_1),\cdots ,{\mathbf{e}}_i(v_n){)}^T$

• 特征向量的第k个元素表示的是波${\mathbf{e}}_i$在顶点$v_k$的振幅；波的频率由对应的特征值${\lambda }_i$决定
• $\mathbf{L}$的特征向量因此被称为自然振动，固有振动（natural vibrations）；特征值是**自然频率，固有频率****（natural frequencies）
• $\mathbf{L}$是对称的，半正定的，它的特征向量构成了一个${\mathbb{R}}^n$正交基，因此，每个向量$\mathbf{f}=(f_1,\cdots ,f_n{)}^T$可以表示为$\mathbf{f}={\sum }_{i=1}^n<{\mathbf{e}}_i,\mathbf{f}>{\mathbf{e}}_i$

利用这个表示可以重建滤波后的函数，例如只使用前m个特征向量$\tilde{\mathbf{f}}={\sum }_{i=1}^m<{\mathbf{e}}_i,\mathbf{f}>{\mathbf{e}}_i$

• 这种方法比较昂贵，因为需要计算特别大的Laplace矩阵的特征向量
• diffusion flow是一种便宜的方法

4.2 扩散流（Diffusion flow）

很多物理过程可以描述为扩散流（例如热扩散 heat diffusion和布朗运动）。扩散流可用扩散方程表示：
$$\frac{\partial f(\mathbf{x},t)}{\partial t}=\lambda \Delta f(\mathbf{x},t)$$
把Laplace算子替换成流形上的Laplace-Beltrami算子，然后对空间和时间离散化：

• 空间离散化，顶点的采样值替代函数，这样，每个顶点得到一个进化方程：

$$\frac{\partial }{\partial t}f(v_i,t)=\lambda \Delta f(v_i,t),i=1,\cdots ,n$$

也可写成矩阵形式：
$$\frac{\partial \mathbf{f}(t)}{\partial t}=\lambda \mathbf{L}\mathbf{f}(t)$$

• 时间离散化，使用时间步长：

$$\mathbf{f}(t+h)=\mathbf{f}(t)+h\frac{\partial \mathbf{f}(t)}{\partial t}=\mathbf{f}(t)+h\lambda \mathbf{L}\mathbf{f}(t)$$

上式也叫显式欧拉积分；也可使用隐式欧拉积分：
$$(\mathbf{I}\mathbf{d}-h\lambda \mathbf{L})\mathbf{f}(t+h)=\mathbf{f}(t)$$

隐式积分计算比显示积分麻烦，但是数值稳定性高

为了平滑网格，简单应用上述规则到每个顶点位置$({\mathbf{x}}_1,\cdots ,{\mathbf{x}}_n{)}^T$:
$${\mathbf{x}}_i\leftarrow {\mathbf{x}}_i+h\lambda \Delta {\mathbf{x}}_i$$
因为顶点位置的Laplace-Beltrami与平均曲法向有关$\Delta \mathbf{x}=-2H\mathbf{n}$，上述方法叫做Laplacian平滑，或平均曲率流

• 沿着法向运动仅对余切拉普拉斯有效；统一拉普拉斯不一定成立。统一拉普拉斯试图把每个顶点移到它的1-邻域的重心上。这种平滑导致切线上的松弛，如下图。这种平滑可能是需要的特色，也可能是一种缺陷。

左图是原始网格，中图是uniform Laplacian平滑的结果，右图是cotangent Laplacian的结果

• 最后，注意到高阶Laplacian流$\partial \mathbf{f}/\partial t=\lambda {\Delta }^k\mathbf{f}$也可以用。其中，高阶拉普拉斯可以使用递归求解${\Delta }^{k}f=\Delta ({\Delta }^{k-1}f)$
• 高阶拉普拉斯效果更好，但是计算开销更大。2阶拉普拉斯能取得$C^1$平滑，拉普拉斯只能达到$C^0$边缘平滑

4.3 光顺

不用的应用对光顺的度量方式不同，不过大体遵循最简形状原则（principle of simplest shape）：表面应该与无关紧要的细节或振动（oscillations）无关

• 可以通过合适的能量来描述

• 常用的光顺函数是膜能量（membrane energy）

$$E_M(\mathbf{x})={\iint }_{\Omega }\sqrt{det(\mathbf{I})}dudv$$

这个能量度量了表面面积，对于修补表面边缘$\partial \Omega$来说，最小面积对应的就是膜表面或称为最小表面

• 然而，这个能量方程高度非线性，还包含了第一基本型，最小化计算困难。使用一阶偏导替代第一基本型，得到狄利克雷能量：

$${\tilde{E}}_M(\mathbf{x})={\iint }_{\Omega }||{\mathbf{x}}_u|{|}^2+||{\mathbf{x}}_v|{|}^{2}dudv$$

为求最小化，使用变分法（calculus of variations）

• 把变分法从1D扩展到流形上，得到流形上的欧拉-拉格朗日方程：

$$\Delta \mathbf{x}(u,v)=0for(u,v)\in \Omega$$

• 离散化：
  • 用顶点坐标组成的向量替代连续坐标
  • 使用离散Laplace-Beltrami算子替代Laplace

$$\mathbf{L}\mathbf{x}=0$$

• 如果目标是最小化曲率而不是表面积，使用thin-plate能量：

$$E_{TP}(\mathbf{x})={\iint }_{\Omega }{\kappa }_1^2+{\kappa }_2^{2}dudv$$

将其线性化为二阶微分量：
$${\tilde{E}}_{TP}(\mathbf{x})={\iint }_{\Omega }||{\mathbf{x}}_{uu}|{|}^2+2||{\mathbf{x}}_{uv}|{|}^2+||{\mathbf{x}}_{vv}|{|}^{2}dudv$$
对应的Euler-Lagrange方程式二阶形式${\Delta }^2\mathbf{x}(u,v)=0$，进一步表示为二阶Laplacian系统${\mathbf{L}}^2\mathbf{x}=0$

• 还可以使用更高阶的PDE，比如最小变分表面 ${\Delta }^3\mathbf{x}=0$

• Uniform Laplacian在顶点密集的地方变化出现误差；contangent能确保正确的结果