Geometry Processing

Acquisition and Processing of 3D Geometry。参考书：Polygon Mesh Processing。下文章节号为课程编号，括号内标注对应原书章节。

Chapter 1: Introduction（原书第 2 章）

Half-edge Data Structure

A half-edge stores a reference to its twin, as well as references to the previous and next half-edges along the same face or hole, plus the origin of the half-edge and its incident face.
A vertex stores its position and a reference to an arbitrary half-edge that originates from that vertex.
A face stores an arbitrary half-edge belonging to that face. A half-edge data structure stores arrays of vertex, face, and half-edge records.

遍历一个面

start_he = f.halfedge;
he = start_he;
do {
  he = he.next;
} while (he != start_he);

遍历一个顶点（逆时针）

start_he = v.halfedge;
he = start_he;
do {
  he = he.prev.twin;
} while (he != start_he);

遍历一个顶点（顺时针）

start_he = v.halfedge;
he = start_he;
do {
  he = he.twin.next;
} while (he != start_he);

Edge Flip

Basic concepts:
  edge.previous       ---->(becomes) edge.next
  edge.twin.next      ---->(becomes) edge.previous

So the same as its twin:
  edge.twin.previous  ---->(becomes) edge.twin.next
  edge.next           ---->(becomes) edge.twin.previous

Dual Graph

顶点变成面，面变成顶点。

Euler–Poincaré Formula

\[ V - E + F = 2(1 - g) \]

对于 genus 为 \(g\) 的闭合多边形网格，顶点数 \(V\)、边数 \(E\)、面数 \(F\) 之间满足上面的欧拉公式。

对于三角网格：

\[ \begin{cases} F \approx 2V \\ E \approx 3V \\ \text{Average valence} = 6 \end{cases} \]

对于四边形网格：

\[ \begin{cases} F \approx V \\ E \approx 2V \\ \text{Average valence} = 4 \end{cases} \]

Chapter 2: Registration

Iterative Closest Points（ICP）算法

1. 选取点的子集 \(\mathbf{p}_i\)

Stable sampling
Normal-based sampling
Slippage analysis

2. 把每个 \(\mathbf{p}_i\) 匹配到另一个 scan 上最近的点 \(\mathbf{q}_i\)

使用层次化的 BSP tree（也可以用 KD-tree）

3. 剔除「坏」的配对 \((\mathbf{p}_i,\mathbf{q}_i)\)

只匹配兼容的点

4. 计算旋转 \(R\) 和平移 \(t\) 使下式最小

\[ \begin{cases} \min_{R,t} \sum\|\mathbf{p}_i - R\mathbf{q}_i - \mathbf{t}\|^2 & \textrm{Point-to-Point} \\[4pt] \min_{R,t} \sum\big[(\mathbf{p}_i - R\mathbf{q}_i - \mathbf{t})^T \mathbf{n}_p\big]^2 & \textrm{Point-to-Plane} \end{cases} \]

约束为 \(R^TR = I,\ \det(R) = 1\)。

Point-to-Point 情形有闭式解：

先计算 \(P\) 和 \(Q\) 的质心：

\[ \begin{cases} \bar{p} = \dfrac{\sum p_i}{n} \\[6pt] \bar{q} = \dfrac{\sum q_i}{n} \end{cases} \]

去中心化：

\[ \begin{cases} \tilde{p}_i = p_i - \bar{p}\\ \tilde{q}_i = q_i - \bar{q} \end{cases} \]

令 \(\tilde{P} = [\tilde{p}_1 \ \ldots \ \tilde{p}_n]\)，\(\tilde{Q}\) 同理。对 \(\tilde{Q}\tilde{P}^T\) 做 SVD：

\[ \tilde{Q}\tilde{P}^T = U\Sigma V^T \]

于是

\[ R = V \begin{bmatrix} 1 & & \\ & 1 & \\ & & \det(VU^T) \end{bmatrix} U^T, \qquad \mathbf{t} = \bar{p} - R\bar{q} \]

5. 对齐后迭代 \(\mathbf{q}_i \leftarrow R\mathbf{q}_i + \mathbf{t}\)

方法一：固定迭代次数
方法二：迭代直到 loss 小于阈值
方法三：两者结合

Chapter 3: Explicit and Implicit（原书第 1 章）

Explicit vs. Implicit

Zero level set：满足 \(F(x) = 0\) 的点集。

Implicit

Signed Distance Function：

\[ F(x,y,z) = \begin{cases} > 0 & \text{outside the surface}\\ = 0 & \text{on the surface}\\ < 0 & \text{inside the surface} \end{cases} \]

平面：

\[ F(\mathbf{p}) = \mathbf{n}^T\mathbf{p} + d \]

Explicit

Curve：\(C(t)\)，一个参数
Surface：\(S(u,v)\)，两个参数

\[ \vec{p} = \lambda_1\vec{a} + \lambda_2\vec{b} + \lambda_3\vec{c} \quad \text{with} \quad \lambda_1 + \lambda_2 + \lambda_3 = 1 \]

Implicit → Explicit：Marching Cubes / Squares

Classify grid nodes as inside/outside
Classify cell：\(2^8\) 种配置
- In/out for each corner（cube 有 8 个顶点）
Compute intersection points
- Linear interpolation along edges
Connect them by triangles
- 每种配置查表（look-up table）
- 通过修改后的表来消除歧义

Explicit → Implicit：Fast Marching

在 mesh 邻域内用精确距离初始化
向外 fast-march
向内 fast-march

Signed Distance Heuristic

Closest point：\(\mathbf{p} = \alpha\mathbf{p}_i + (1 - \alpha)\mathbf{p}_j\)
Interpolated normal：\(\mathbf{n} = \alpha \mathbf{n}_i + (1 - \alpha)\mathbf{n}_j\)
Inside if \((\mathbf{q} - \mathbf{p})^T\mathbf{n} < 0\)

Chapter 4 & 5: Differential Geometry（原书第 3 章）

点云的法向估计

1. Least-squares estimation

平面定义为 \(ax + by + cz + d = 0\)，平面的法向是 \((a,b,c)\)。\(d\) 可以取任意非零值，因为它不影响法向的方向（通常取 \(-1\)，于是方程变成 \(A\mathbf{n} = \vec{1}\)）。

2. PCA based normal estimation

对于三维点云，法向是协方差矩阵最小特征值对应的特征向量。步骤：

计算点云的均值（按维度）
减去均值
计算协方差矩阵 \(C = X^TX\)（未做 \(1/n\) 归一化，不影响特征向量方向）
对协方差矩阵做特征分解，或对 \(X\) 做 SVD
取特征值按降序排列后的第三个特征向量（2D 情形取第二个），即最小特征值对应的方向

这里 \(W\) 是 \(C\) 的特征向量矩阵：因为对 \(X\) 做 SVD 时，\(U\) 的列是 \(XX^T\) 的特征向量，\(W\) 的列是 \(X^TX\) 的特征向量。

\[ C = \begin{cases} Q\Lambda Q^T & \textrm{Eigen Decomposition} \\ W\Sigma^2 W^T & \textrm{SVD Decomposition} \end{cases} \qquad \text{其中}\ X = U\Sigma W^T \]

为什么需要微分几何

需要计算：

Surface curvature
Parameterization distortion
Deformation energies

微分几何的分类：

Local：研究曲线和曲面的局部性质，只取决于某点邻域内的行为。
Global：研究整体性质。

Parametric Curves

定义：

\[ \mathbf{x}(t) = \begin{pmatrix} x(t) \\ y(t) \\ z(t) \end{pmatrix} \]

于是

\[ \mathbf{x}_t(t) := \frac{d\mathbf{x}(t)}{dt} = \begin{pmatrix} dx(t)/dt\\ dy(t)/dt\\ dz(t)/dt \end{pmatrix} \]

令 \(t_i = a + i\Delta t\)，曲线为 \(\alpha(t_i)\)，则弧长为

\[ L(c) = \lim_{\Delta \rightarrow 0}\sum_{t_i}\frac{\|\alpha(t_i) - \alpha(t_i + \Delta)\|}{\Delta}\Delta = \int_a^b \|\alpha^{\prime}(t)\|\,dt \]

对于弧长参数化，有

\[ \|\alpha^{\prime}(s)\| = 1 \]

它把参数区间 \([a,b]\) 映射到 \([0,L]\)，其中 \(L = l(a,b) = \int^b_a\|\mathbf{x}^{\prime}(u)\|\,du\)。

The Frenet Frame

由 \(\|\alpha^{\prime}(s)\| = 1\) 得 \(\|\alpha^{\prime}(s)\|^2 = 1\)，两边对 \(s\) 求导得

\[ 2\,\alpha^{\prime}(s)\cdot\alpha^{\prime\prime}(s) = 0 \]

\(\alpha^{\prime}(s)\) 是 \(s\) 处的切方向，记作 \(\mathbf{t}\)。
\(\kappa(s) = \|\alpha^{\prime\prime}(s)\|\) 是 \(\alpha\) 在 \(s\) 处的曲率。
- 直观上，曲率度量一条曲线偏离直线的程度。
- 由 \(\alpha^{\prime}(s)\cdot\alpha^{\prime\prime}(s) = 0\) 可知，二阶导数与曲线的一阶导数（切方向）正交，因此二阶导数在法方向上，于是 \(\alpha^{\prime\prime}(s) = \kappa(s)\,\mathbf{n}(s)\)。
- 曲率半径定义为曲率的倒数：\(R(s) = 1/\kappa(s)\)。
\(\mathbf{b}(s) = \mathbf{t}(s) \times \mathbf{n}(s)\) 是副法方向（binormal）。
\(\mathbf{t},\mathbf{n},\mathbf{b}\) 一起构成 Frenet Frame。
\(\mathbf{t}(s),\mathbf{n}(s)\) 张成 osculating plane（密切平面）。
\(\mathbf{n}(s),\mathbf{b}(s)\) 张成 normal plane（法平面）。
\(\mathbf{t}(s),\mathbf{b}(s)\) 张成 rectifying plane（从切平面）。
满足 \(\alpha^{\prime}(s) = 0\) 的点称为 0 阶奇点。
满足 \(\alpha^{\prime\prime}(s) = 0\) 的点称为 1 阶奇点。

三个向量之间的关系（Frenet–Serret 公式）：

\[ \begin{cases} \mathbf{t}^{\prime} = \kappa\,\mathbf{n} \\ \mathbf{n}^{\prime} = -\kappa\,\mathbf{t} + \tau\,\mathbf{b} \\ \mathbf{b}^{\prime} = -\tau\,\mathbf{n} \end{cases} \]

其中 \(\mathbf{b}^{\prime}\) 度量了曲线在 \(s\) 处脱离密切平面的快慢；\(\tau\) 是 \(\alpha\) 在 \(s\) 处的扭率（torsion）。

Curvature：偏离直线的程度
Torsion：偏离平面的程度

小结

求什么	怎么求
\(\mathbf{t}(s)\)	\(\mathbf{t} = \alpha^{\prime} = (x',y',z')\)
\(\mathbf{n}(s)\)	由 \(\mathbf{t}' = \kappa\mathbf{n}\) 得到
\(\mathbf{b}(s)\)	\(\mathbf{b} = \mathbf{t} \times \mathbf{n}\)
\(\kappa(s)\)	\(\kappa = \\|\mathbf{t}'\\|\)
\(\tau(s)\)	由 \(\mathbf{b}' = -\tau\mathbf{n}\) 得到

Cauchy–Crofton Theorem

\[ L(c) = \frac12\iint n(p,\theta)\,dp\,d\theta \cong \frac{rn\pi}{8} \]

其中 \(n\) 是相交次数。构造一族间距为 \(r\)、平行于 \(x\) 轴的直线，把这族直线绕原点分别旋转 \(\frac{\pi}{4},\frac{\pi}{2},\frac{3\pi}{4}\)，就得到四族平行线。

Parametric Surfaces

\[ \mathbf{x}(u,v) = \begin{pmatrix}x(u,v)\\y(u,v)\\z(u,v)\end{pmatrix} \]

uv 平面上的曲线 \([u(t),v(t)]\) 定义了曲面 \(\mathbf{x}(u,v)\) 上的一条曲线：

\[ \mathbf{c}(t) = \mathbf{x}(u(t), v(t)) \]

法向量：

\[ \mathbf{n} = \frac{\mathbf{x}_u \times \mathbf{x}_v}{\|\mathbf{x}_u \times \mathbf{x}_v\|} \]

因为 \(\|\mathbf{n}(u,v)\| = 1\)，所以 \(\mathbf{n}(u,v) \cdot \mathbf{n}'(u,v) = 0\)，即法向的变化位于切平面内。

注意 \(\mathbf{x}_u, \mathbf{x}_v\) 不一定正交，但要求 \(\mathbf{x}_u \times \mathbf{x}_v \neq 0\)：

\[ \frac{\partial\mathbf{x}(u,v)}{\partial u} = \mathbf{x}_u = \begin{pmatrix} \partial x/\partial u \\ \partial y/\partial u \\ \partial z/\partial u \end{pmatrix}, \qquad \frac{\partial\mathbf{x}(u,v)}{\partial v} = \mathbf{x}_v = \begin{pmatrix} \partial x/\partial v \\ \partial y/\partial v \\ \partial z/\partial v \end{pmatrix} \]

所以 \(\mathbf{J} = [\mathbf{x}_u,\mathbf{x}_v]\) 是一个 \(3 \times 2\) 矩阵，这两列定义了切平面：

\[ \frac{\mathrm{d}\alpha(t)}{\mathrm{d}t} =\frac{\mathrm{d}\mathbf{x}(u(t),v(t))}{\mathrm{d}t} =\frac{\partial\mathbf{x}}{\partial u}\frac{\mathrm{d}u}{\mathrm{d}t} +\frac{\partial\mathbf{x}}{\partial v}\frac{\mathrm{d}v}{\mathrm{d}t} =\mathbf{x}_u u_t+\mathbf{x}_v v_t \]

First Fundamental Form

\[ \mathbf{I}=\mathbf{J}^T\mathbf{J} =\begin{bmatrix}E&F\\F&G\end{bmatrix} :=\begin{bmatrix}\mathbf{x}_u^T\mathbf{x}_u&\mathbf{x}_u^T\mathbf{x}_v\\\mathbf{x}_u^T\mathbf{x}_v&\mathbf{x}_v^T\mathbf{x}_v\end{bmatrix} \]

由它可以得到：

Angles：\(\cos\theta = \dfrac{F}{\sqrt{EG}}\)
Length：\(l(a,b) = \displaystyle\int_a^b\sqrt{Eu_t^2+2Fu_tv_t+Gv_t^2}\,\mathrm{d}t\)
Area：\(A = \displaystyle\iint_U\sqrt{EG-F^2}\,\mathrm{d}u\,\mathrm{d}v\)

Second Fundamental Form

\[ \mathbf{II}=\begin{pmatrix}e&f\\f&g\end{pmatrix} :=\begin{pmatrix}\mathbf{x}_{uu}^T\mathbf{n}&\mathbf{x}_{uv}^T\mathbf{n}\\\mathbf{x}_{uv}^T\mathbf{n}&\mathbf{x}_{vv}^T\mathbf{n}\end{pmatrix} \]

Normal curvature \(\kappa_n(t)\) 定义为法截线 \(\mathbf{c}(t)\) 在点 \(\mathbf{p} = \mathbf{x}(u,v)\) 处的曲率，可以由两个基本形式算出：

\[ \kappa_n(\bar{\mathbf{t}}) =\frac{\bar{\mathbf{t}}^T\mathbf{II}\,\bar{\mathbf{t}}}{\bar{\mathbf{t}}^T\mathbf{I}\,\bar{\mathbf{t}}} =\frac{eu_t^2+2fu_tv_t+gv_t^2}{Eu_t^2+2Fu_tv_t+Gv_t^2} \]

其中 \(\mathbf{t}=u_t\mathbf{x}_u+v_t\mathbf{x}_v\)，\(\bar{\mathbf{t}} = (u_t,v_t)\)，\(\mathbf{t}\) 是曲面点 \(\mathbf{p}\) 处的切向量。

曲面的曲率性质可以通过考察 \(\mathbf{p}\) 处所有法截线的曲率来刻画，也就是让切向量 \(\mathbf{t}\) 绕曲面法向旋转一圈。

Principal curvatures

旋转过程中会出现两个不同的极值：

Maximum curvature：\(\kappa_1 = \max_{\psi}\kappa_n(\psi)\)
Minimum curvature：\(\kappa_2 = \min_{\psi}\kappa_n(\psi)\)
Euler theorem：\(\kappa_n(\bar{\mathbf{t}}) = \kappa_1\cos^2\psi + \kappa_2\sin^2\psi\)，其中 \(\psi\) 是 \(\mathbf{t}\) 与 \(\mathbf{t}_1\) 的夹角。这说明曲面的曲率完全由两个主曲率决定。
对应的主方向 \(\mathbf{e}_1,\mathbf{e}_2\) 相互正交（也可由 Euler theorem 推出）。
主曲率是形状算子 \(S = \mathbf{I}^{-1}\mathbf{II}\) 的特征值，主方向是它的特征向量。

Special curvatures

Mean curvature：\(H = \dfrac{\kappa_1 + \kappa_2}{2}\)，主曲率的平均
Gaussian curvature：\(K = \kappa_1 \cdot \kappa_2\)，主曲率的乘积

用形状算子表示：

\[ H = \frac12\operatorname{trace}\big(\mathbf{I}^{-1}\mathbf{II}\big), \qquad K = \det\big(\mathbf{I}^{-1}\mathbf{II}\big) = \frac{\det \mathbf{II}}{\det \mathbf{I}} = \frac{eg-f^2}{EG-F^2} \]

Gaussian curvature 可以把曲面上的点分成几类：

\[ \begin{array}{ll} \text{Elliptic} & \text{if } K > 0\\ \text{Hyperbolic} & \text{if } K < 0\\ \text{Parabolic} & \text{if } K = 0\\ \text{Umbilic} & \text{if } \kappa_1 = \kappa_2 \end{array} \]

若处处 \(H = 0\) \(\rightarrow\) minimal surface
若处处 \(K = 0\) \(\rightarrow\) developable surface

Intrinsic Geometry

只依赖第一基本形式的曲面性质：

Length
Angle
Gaussian curvature（Theorema Egregium）——曲面的高斯曲率完全由长度和角度决定

\[ K = \lim_{r \rightarrow 0} \frac{6\pi r - 3C(r)}{\pi r^3} \]

这说明高斯曲率在局部等距变换下不变，因此也是曲面的内蕴量。这里 \(r\) 是曲面上某点周围小邻域的半径，\(C(r)\) 是该点处半径为 \(r\) 的测地圆的周长。

曲面上一点的分类

Isotropic：\(\kappa_1 = \kappa_2\)，各方向相同
Anisotropic：\(\kappa_1 \neq \kappa_2\)，存在不同的主方向

Gauss–Bonnet Theorem

对任意闭合流形曲面，其欧拉特征数 \(\chi = 2 - 2g\)，则

\[ \int_\Omega K(u,v)\,\mathrm{d}u\,\mathrm{d}v = 2\pi\chi \]

Differential Operators

Gradient

\[ \nabla f:=\left(\frac{\partial f}{\partial x_1},\ldots,\frac{\partial f}{\partial x_n}\right) \]

Divergence

\[ \operatorname{div}F=\nabla\cdot F:=\frac{\partial F_1}{\partial x_1}+\ldots+\frac{\partial F_n}{\partial x_n} \]

Chapter 6: Smoothing（原书第 4 章）

1. Differential Operators

\[ \operatorname{div} g:=\left\langle\left(\frac{\partial}{\partial x_1},\ldots\right),\,g\right\rangle=\nabla\cdot g \]

\[ \operatorname{div}\mathbf{G} :=\left\langle\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}\right),\,(G_x,G_y,G_z)\right\rangle =\frac{\partial G_x}{\partial x}+\frac{\partial G_y}{\partial y}+\frac{\partial G_z}{\partial z} \]

与 gradient 得到向量不同，divergence 把各个分量加起来得到标量。

2. Laplace Operator

\[ \operatorname{div}\nabla f=\sum_i\frac{\partial^2f}{\partial x_i^2}=\Delta f \]

其中 \(\Delta = \nabla \cdot \nabla\) 称为 Laplace 算子（梯度的散度），\(f\) 是欧氏空间中的标量函数。

3. Laplace–Beltrami Operator

它把 Laplace 算子推广到定义在曲面上的函数：

\[ \Delta_{\mathcal{S}}\mathbf{x}=\operatorname{div}_{\mathcal{S}}\nabla_{\mathcal{S}}\mathbf{x} = -2H\mathbf{n} \]

其中 \(\Delta_{\mathcal{S}}\) 称为 Laplace–Beltrami 算子，\(H\) 是平均曲率，\(\mathbf{n}\) 是曲面法向。这个关系式在后面反复用到，下文称之为「\(\Delta_{\mathcal{S}}\mathbf{x} = -2H\mathbf{n}\) 关系」。

4. Discrete Curvature

由于多边形网格是分片线性曲面，无法直接套用上述连续概念。因此以下离散微分算子的定义，都基于「网格可解释为光滑曲面的分片线性近似」这一假设。

一般思路是把离散微分性质计算为网格上某点 \(x\) 的局部邻域 \(N(x)\) 上的空间平均。

顶点 one-ring 邻域上的三种平均区域：

Barycentric Cells：连接三角形重心和边的中点

Voronoi Cells：连接外心（三边垂直平分线的交点）

Mixed Cells：对于外心落在 Voronoi cell 之外的情形，用中心顶点对边的中点来替代

5. Discrete Laplace–Beltrami（per vertex）

Uniform Discretization

\[ \Delta f(v_i)=\frac{1}{|\mathcal{N}_1(v_i)|}\sum_{v_j\in\mathcal{N}_1(v_i)}\big(f(v_j)-f(v_i)\big) \]

求和遍历所有 one-ring 邻居 \(v_j \in \mathcal{N}_1(v_i)\)，\(|\mathcal{N}_1(v_i)|\) 是 valence，下标 1 表示 1-ring。等价地也可以写成

\[ \Delta f(v_i)=\frac{1}{|\mathcal{N}_1(v_i)|}\sum_{v_j\in\mathcal{N}_1(v_i)}f(v_j) \;-\; f(v_i) \]

图中蓝线表示 \(\Delta_{\mathcal{S}} \mathbf{x}\) 或 \(\Delta_{\text{uni}}\mathbf{x}\)。

然而，即使顶点是共面配置，得到的向量也可能非零；而根据 \(\Delta_{\mathcal{S}}\mathbf{x} = -2H\mathbf{n}\) 关系，此时它应该为零。这说明 uniform Laplacian 对非均匀网格不是一个合适的离散化。

Cotangent Discretization

\[ \int_{A_i}\operatorname{div}\mathbf{F}(\mathbf{u})\,\mathrm{d}A=\int_{\partial A_i}\mathbf{F}(\mathbf{u})\cdot\mathbf{n}(\mathbf{u})\,\mathrm{d}s \]

\(A_i\) 是平均面积。这是一个更精确的离散化，上式可以化简为

\[ \Delta f(v_i) := \frac{1}{2A_i}\sum_{v_j\in\mathcal{N}_1(v_i)}\big(\cot\alpha_{i,j}+\cot\beta_{i,j}\big)\big(f_j-f_i\big) \]

与 uniform 版本相比，它多了权重 \((\cot\alpha_{i,j} +\cot\beta_{i,j})\)。

6. Discrete Curvature

由 \(\Delta_{\mathcal{S}} \mathbf{x} = -2H\mathbf{n}\) 关系可得：

Mean Curvature（绝对值）

\[ H = \frac{1}{2}\|\Delta_{\mathcal{S}} \mathbf{x}\| \]

Gaussian curvature

\[ K(v_i)=\frac{1}{A_i}\left(2\pi-\sum_{v_j\in\mathcal{N}_1(v_i)}\theta_j\right) \]

有了这两个量，主曲率可以算出来：

\[ \kappa_{1,2}(v_i)=H(v_i)\pm\sqrt{H(v_i)^2-K(v_i)} \]

7. Mesh Quality Criteria

Smoothness
Fairness
Adaptive tessellation（镶嵌）
Triangle shape

8. Laplace Operator on Surfaces

\[ \begin{aligned} &\Delta = M^{-1}C\\ &M = M^T, \quad C = C^T \\ &\Delta^T \neq \Delta \\ &\text{其中}\\ &c_{ij}=\frac12\big(\cot(\alpha_{ij})+\cot(\beta_{ij})\big) \\ &c_{ii}=-\sum_{j\in \mathcal{N}_1(i)}c_{ij} \\ &\sum_j c_{ij}=0 \end{aligned} \]

9. Common Problem Types

Laplace equation（Poisson 方程的特例）

\[ \begin{aligned} & \Delta f = 0 \\ & M^{-1}Cf = 0 \Rightarrow Cf = 0 \end{aligned} \]

边界条件为 \(f(x) = f_0,\ x \in \partial B\)。

Poisson equation

\[ \begin{aligned} & \Delta f = g \\ & M^{-1}Cf = g \Rightarrow Cf = Mg \end{aligned} \]

这里 \(C\) 是对称的，可以把它看成 \(Ax = b\)，其中 \(b = Mg\)。

Diffusion equation

\[ f_t = \frac{\partial f}{\partial t} = \Delta f \]

物理含义：函数 \(f\) 在某一点上随时间的变化率（例如温度或浓度随时间的变化），与该点上函数的空间曲率（例如温度或浓度空间分布的弯曲程度）成正比。

Eigen / spectral analysis

\[ \begin{aligned} & \Delta \phi_i = \lambda_i\phi_i \\ & M^{-\frac{1}{2}}CM^{-\frac{1}{2}}M^{\frac{1}{2}}\phi_i=\lambda_i M^{\frac{1}{2}}\phi_i \\ & \left(M^{-\frac12}CM^{-\frac12}\right)\xi_i=\lambda_i\xi_i \quad \Rightarrow \quad \phi_i = M^{-\frac12}\xi_i \end{aligned} \]

其中 \(\left(M^{-\frac12}CM^{-\frac12}\right)\xi_i=\lambda_i\xi_i\) 说明 \(\xi_i\) 是矩阵 \(M^{-\frac12}CM^{-\frac12}\) 的特征向量，\(\lambda_i\) 是对应的特征值。

10. Spectral Analysis

Fourier Transform：

\[ F(\omega)=\int_{-\infty}^{\infty}f(x)\,\mathrm{e}^{-2\pi\mathrm{i}\omega x}\,\mathrm{d}x = \big\langle f(x),\, \mathrm{e}^{2\pi \mathrm{i}\omega x}\big\rangle \]

其中内积定义为

\[ \langle f,g\rangle = \int_{-\infty}^{\infty} f(x)\overline{g(x)}\,\mathrm{d}x \]

横线代表复共轭。

类比线性代数：

\[ \begin{aligned} &Ax = b \\ &Ae_i = \lambda_i e_i \quad \text{where } \|e_i\| = 1 \\ &x = \sum a_i e_i = \sum\langle e_i,x\rangle e_i \end{aligned} \]

而逆变换

\[ f(x)=\int_{-\infty}^{\infty}F(\omega)\,\mathrm{e}^{2\pi\mathrm{i}\omega x}\,\mathrm{d}\omega =\int_{-\infty}^{\infty}\big\langle f,\mathrm{e}^{2\pi\mathrm{i}\omega x}\big\rangle\,\mathrm{e}^{2\pi\mathrm{i}\omega x}\,\mathrm{d}\omega \]

和上面的展开式是同一件事：都是把函数在一组正交基上展开再重组。

Fourier Analysis on Meshes：

\[ \Delta\big(\mathrm{e}^{2\pi\mathrm{i}\omega x}\big) = \frac{\mathrm{d}^2}{\mathrm{d}x^2}\mathrm{e}^{2\pi\mathrm{i}\omega x}=-\left(2\pi\omega\right)^2\mathrm{e}^{2\pi\mathrm{i}\omega x} \]

形式上类似 \(Ax = \lambda x\)，所以基函数正是 Laplacian 算子的特征向量。

由于我们知道如何在离散三角网格上离散化 Laplace–Beltrami 算子，就可以利用这一点来定义三角网格上的傅里叶变换：找到网格上函数空间的一组基，这组基是离散 Laplace–Beltrami 算子的特征向量。这些特征向量可以用来分析和重构定义在三角网格上的函数，与傅里叶变换在欧氏空间中的应用类似。

11. Discrete Laplace–Beltrami（per mesh）

在三角网格上离散化时，把连续函数 \(f(x)\) 替换为在 \(n\) 个网格顶点上的采样值向量：

\[ f:\mathcal{S}\to\mathbb{R}\quad\longrightarrow\quad\big(f(v_1),\ldots,f(v_n)\big)^T \]

离散 Laplace–Beltrami 的稀疏矩阵形式：

\[ \begin{pmatrix}\Delta f(v_1)\\\vdots\\\Delta f(v_n)\end{pmatrix} =\mathbf{L}\begin{pmatrix}f(v_1)\\\vdots\\f(v_n)\end{pmatrix} \]

其中 \(\mathbf{L} = M^{-1}C \in \mathbb{R}^{n\times n}\)：

\[ \begin{aligned} &\mathbf{C}_{ij}= \begin{cases} \cot\alpha_{ij}+\cot\beta_{ij}, & i\neq j,\ j\in\mathcal{N}_1(v_i)\\ -\sum_{v_j\in\mathcal{N}_1(v_i)}(\cot\alpha_{ij}+\cot\beta_{ij}) & i=j\\ 0 & \mathrm{otherwise} \end{cases}\\[6pt] &\mathbf{M}^{-1}=\operatorname{diag}\left(\ldots,\frac{1}{2A_i},\ldots\right) \end{aligned} \]

连续情形下 Laplace–Beltrami 算子的特征函数 \(e_{\omega}(x)\)，现在变成 Laplace 矩阵的特征向量 \(\mathbf{e}_1,\ldots,\mathbf{e}_n\)：一个 \(n\) 维特征向量 \(\mathbf{e}_i\) 可以看作连续特征函数 \(\mathbf{e}_i(\mathbf{x})\) 的离散采样 \((\mathbf{e}_i(v_1),\ldots,\mathbf{e}_i(v_n))^T\)。\(\mathbf{e}_i\) 的第 \(k\) 个分量对应波 \(\mathbf{e}_i\) 在顶点 \(v_k\) 处的振幅，波的频率由对应的特征值 \(\lambda_i\) 决定。因此 \(\mathbf{L}\) 的特征向量被称为「natural vibrations」，特征值被称为「natural frequencies」。

12. Spectral Analysis 流程

构造 Laplace–Beltrami 矩阵 \(\mathbf{L} = \Delta\)
计算最小的 \(k\) 个特征值所对应的特征向量 \(\mathbf{e}_1,\ldots,\mathbf{e}_k\)
用这些特征向量按分量重构网格

13. Diffusion Flow

它相当于用一个高斯核去衰减高频，而不是像 spectral analysis 那样把阈值 \(\omega_{\max}\) 以上的频率直接截断。

以下全部基于扩散方程。

1. Diffusion Flow on Height Fields

\[ \frac{\partial f}{\partial t}=\lambda\Delta f \]

其中 \(\lambda\) 是扩散系数。

2. Diffusion Flow on Vertex（空间离散化）

\[ \frac{\partial}{\partial t}f(v_i,t) = \lambda\Delta f(v_i,t),\quad i=1,\ldots,n \]

3. Diffusion Flow on Mesh（空间离散化）

矩阵形式：

\[ \frac{\partial \mathbf{f}(t)}{\partial t} = \lambda \mathbf{L}\mathbf{f}(t) = \lambda \Delta_{\mathcal{S}} \mathbf{f} = -2\lambda H \mathbf{n} \]

4. Explicit Euler Integration（时间离散化）

\[ \begin{aligned} \mathbf{p}_i^{(t+1)}&=\mathbf{p}_i^{(t)}+\lambda\Delta\mathbf{p}_i^{(t)}\\ \mathbf{P}^{(t+1)}&=\left(\mathbf{I}+\lambda\mathbf{L}\right)\mathbf{P}^{(t)} \end{aligned} \]

这要求 \(\lambda\) 足够小才能稳定，否则会出现震荡。

5. Implicit Euler Integration（时间离散化）

\[ \begin{aligned} \mathbf{p}_i^{(t+1)}&=\mathbf{p}_i^{(t)}+\lambda\Delta\mathbf{p}_i^{(t+1)} \\ \mathbf{P}^{(t)}&=\left(\mathbf{I}-\lambda\mathbf{L}\right)\mathbf{P}^{(t+1)} \end{aligned} \]

但由于 \(\mathbf{I} - \lambda \mathbf{L}\) 不对称，求解很慢。

解决办法是把 \(\mathbf{L}\) 替换为 \(\mathbf{M}^{-1}\mathbf{C}\)，于是上式变成

\[ (\mathbf{M}-\lambda\mathbf{C})\mathbf{P}^{(t+1)}=\mathbf{M}\mathbf{P}^{(t)} \]

这是一个稀疏对称正定系统，可以用：

迭代共轭梯度法
Sparse Cholesky

14. Energy Minimization

Fairness

Idea：惩罚「不美观的行为」
曲面 fairing 的目标是计算尽可能光滑的形状

1. 度量 fairness

Principle of the simplest shape
Physical interpretation

2. 最小化某个 fairness functional

Surface area、curvature 等。

Membrane energy

\[ \int_{\mathcal{S}} \mathrm{d}A \rightarrow \min \quad \text{with}\ \ \delta \mathcal{S} = \mathbf{c} \]

它度量曲面 \(\mathcal{S}\) 的面积。但 membrane energy 高度非线性，因此常用其线性化形式 Dirichlet energy：

\[ \tilde{E}_{\mathrm{M}}(\mathbf{x})=\iint_{\Omega}\left\|\mathbf{x}_{u}\right\|^{2}+\left\|\mathbf{x}_{v}\right\|^{2}\,\mathrm{d}u\,\mathrm{d}v \]

Thin-plate energy

目标是最小化曲率：

\[ \int_{\mathcal{S}}\kappa_1^2+\kappa_2^2\,\mathrm{d}A \to \min \quad\mathrm{with}\quad \delta\mathcal{S}=\mathbf{c},\quad\mathbf{n}(\delta\mathcal{S})=\mathbf{d} \]

它的线性化形式为

\[ \tilde{E}_{\mathrm{TP}}(\mathbf{x})=\iint_{\Omega}\left\|\mathbf{x}_{uu}\right\|^2+2\left\|\mathbf{x}_{uv}\right\|^2+\left\|\mathbf{x}_{vv}\right\|^2\,\mathrm{d}u\,\mathrm{d}v \]

Calculus of Variations

1D membrane（线性形式）能量：

\[ E(f)=\int_{a}^{b}f'(x)^{2}\,\mathrm{d}x \rightarrow \min \]

约束是固定 \(f(a)\) 和 \(f(b)\) 的边界条件。取到最小值的必要条件是

\[ f'' = \Delta f = 0 \]

这称为 Euler–Lagrange 方程：它表明在最小值处，\(E(f)\) 关于 \(f\) 的一阶变分必须为零。

Bivariate Variational Calculus

Membrane surfaces

把上面 Dirichlet energy 的结论推广过来：

\[ \tilde{E}_\mathrm{M}(\mathbf{x})\to\min \quad\Leftrightarrow\quad \Delta\mathbf{x}(u,v)=0\ \mathrm{for}\ (u,v)\in\Omega \]

再把这个连续形式转到离散三角网格上：

\[ \begin{cases} \mathbf{x}(u,v) \rightarrow \mathbf{x}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\right)^{T} \\ \text{使用离散 Laplace–Beltrami 算子} \end{cases} \]

得到

\[ \mathbf{L}\mathbf{x} = \mathbf{0} \qquad\text{或}\qquad \Delta_{\mathcal{S}}\mathbf{x} = \mathbf{0} \]

Thin-plate surfaces

\[ \Delta_{\mathcal{S}}^2\mathbf{x} = \mathbf{0} \qquad\text{或}\qquad \mathbf{L}^2\mathbf{x} = \mathbf{0} \]

Higher order

最小化曲率的变化：

\[ \iint_{\Omega}\left(\frac{\partial\kappa_1}{\partial\mathbf{t}_1}\right)^2+\left(\frac{\partial\kappa_2}{\partial\mathbf{t}_2}\right)^2\,\mathrm{d}u\,\mathrm{d}v \]

得到

\[ \mathbf{L}^k\mathbf{x}=\mathbf{0} \]

等于 0 说明 fair surface 确实已经 as smooth as possible，不然总可以继续进行 \(k\) 阶 diffusion flow。

Laplacian flow 的稳定曲面：

\[ \frac{\partial\mathbf{x}}{\partial t}=\Delta_{\mathcal{S}}^k\mathbf{x} \]

Appendix

File Type

OFF 格式

第一行（在 vscode 中不显示）：定义图中显示的顶点大小
第二行（vscode 中的第一行）：OFF，固定模板
第三行：定义顶点数、面数、边数（边数通常设为 0，没有影响）
第四行 ~ 第 \((n_{\text{vertices}} + 4)\) 行：顶点表，定义顶点坐标
第 \((n_{\text{vertices}} + 4)\) 行 ~ 第 \((n_{\text{vertices}} + 4 + n_{\text{faces}})\) 行：面表
- 先给出该面包含的顶点数
- 再列出这些顶点的索引