5688 字

28 分钟

计算机图形学笔记(三)：几何建模与处理

2025-07-23

计算机图形学

C++

/

编程

/

计算机图形学

/

OpenGL

NOTE
这一部分内容相对抽象，但贝塞尔曲线在 UI、字体、动画里随处可见——学完之后就明白为什么 Photoshop 钢笔工具那么直觉。重点在贝塞尔 + B 样条 + 网格半边结构这三块，曲率/切向量这些微分几何量在 Part 1 §4.1 已经推过，这里只做必要回顾并交叉引用。

目录#

参数曲线基础 — 参数化优势 / 常见曲线族（与 Part 1 §4.1 互补）
贝塞尔曲线 — 伯恩斯坦基、de Casteljau 算法及其正确性证明（canonical 位置）
B 样条与 NURBS — Cox-de Boor 递推、有理扩展
网格处理 — 半边结构、QEM 简化、Loop 细分、拉普拉斯平滑

十二、参数曲线基础#

TIP
本章与 Part 1 §4.1 的边界：切向量 / 曲率 / 弧长这些微分几何定义已在 Part 1 推过。本章只快速复用，重点放在”为什么参数曲线在图形学里好用”以及常见曲线族的参数形式。

12.1 为什么选参数形式#

12.1.1 三种曲线表示法的对比#

表示法	形式	优点	缺点
显式	$y = f(x)$	求值直接	无法表示垂线/闭合曲线；多值函数不适用
隐式	$F(x,y,z) = 0$	判断”点是否在曲线上”很容易	求曲线上的点要解方程；不自然支持参数化采样
参数	$\mathbf{C}(t) = (x(t), y(t), z(t))$	采样、动画、微分几何计算统一；支持任意拓扑	一个曲线可以有无数种参数化

图形学几乎全用参数形式：动画里 $t$ 直接当时间轴，光栅化时均匀或自适应地在 $[0,1]$ 上采样，导数直接拿来做切向量/速度。

12.1.2 常见曲线的参数方程（速查）#

直线（两点线性插值）：

\mathbf{L}(t) = (1-t)\mathbf{P}_0 + t\mathbf{P}_1, \quad t \in [0, 1]

切向量： $\mathbf{L}'(t) = \mathbf{P}_1 - \mathbf{P}_0$ （常向量）
曲率： $\kappa \equiv 0$

圆（参数化成三角函数）：

\mathbf{C}(t) = \mathbf{c} + r\begin{pmatrix} \cos t \\ \sin t \end{pmatrix}, \quad t \in [0, 2\pi]

$\|\mathbf{C}'(t)\| = r$ ， $\kappa = 1/r$
三维中用两个正交单位向量 $\mathbf{u}, \mathbf{v}$ 张成平面即可推广

椭圆：

\mathbf{E}(t) = \mathbf{c} + (a\cos t, b\sin t), \quad \kappa(t) = \frac{ab}{(a^2\sin^2 t + b^2\sin^2 t)^{3/2}}

长轴端点曲率 $b^2/a^3$ ，短轴端点 $a^2/b^3$ ——长轴端点”更尖”

螺旋线（3D）： $\mathbf{S}(t) = (a\cos t, a\sin t, bt)$

NOTE
曲率公式回顾（Part 1 §4.1.2 已推）： $\kappa(t) = \|\mathbf{C}'(t) \times \mathbf{C}''(t)\| / \|\mathbf{C}'(t)\|^3$ 。平面曲线可简化为 $\kappa = |x'y'' - y'x''| / (x'^2 + y'^2)^{3/2}$ 。

12.1.3 弧长参数化#

弧长函数：

s(t) = \int_a^t \|\mathbf{C}'(\tau)\| \, d\tau

以 $s$ 为参数时 $\|\mathbf{C}'(s)\| = 1$ ，于是 $\kappa(s) = \|\mathbf{C}''(s)\|$ ，几何公式变得非常干净。但对多数曲线（贝塞尔/样条）， $s(t)$ 无解析反函数，实际工程里要靠数值查表实现”按弧长匀速走”。

十三、贝塞尔曲线（canonical）#

WARNING
本章是全系列贝塞尔曲线的规范位置。Part 1 §4.1 和 Part 6 里再提到贝塞尔时，都回链这里。

13.1 伯恩斯坦基函数#

13.1.1 定义#

$n$ 次伯恩斯坦基函数：

B_i^n(t) = \binom{n}{i} t^i (1-t)^{n-i}, \quad i = 0, 1, \ldots, n

其中二项式系数 $\binom{n}{i} = n! / [i!(n-i)!]$ 。

13.1.2 五条核心性质（带证明）#

性质 1 — 非负性： $B_i^n(t) \geq 0$ 在 $t \in [0,1]$ 上。

证：三个因子 $\binom{n}{i}, t^i, (1-t)^{n-i}$ 在该区间上均非负。

性质 2 — 单位分解： $\sum_{i=0}^{n} B_i^n(t) = 1$ 。

证：二项式定理直接给出 $\sum_i \binom{n}{i} t^i (1-t)^{n-i} = (t + (1-t))^n = 1$ 。

几何意义：性质 1 + 2 说明 $\{B_i^n\}$ 是一组凸权重，这是贝塞尔曲线凸包性的根本。

性质 3 — 对称性： $B_i^n(t) = B_{n-i}^n(1-t)$ 。

证： $B_{n-i}^n(1-t) = \binom{n}{n-i}(1-t)^{n-i}t^i = \binom{n}{i}t^i(1-t)^{n-i} = B_i^n(t)$ 。

性质 4 — 端点值： $B_i^n(0) = \delta_{i0}$ ， $B_i^n(1) = \delta_{in}$ 。

这直接给出贝塞尔曲线必过首尾控制点。

性质 5 — 递推关系（de Casteljau 的数学基础）：

B_i^n(t) = (1-t)\, B_i^{n-1}(t) + t\, B_{i-1}^{n-1}(t)

证： $\binom{n}{i} = \binom{n-1}{i} + \binom{n-1}{i-1}$ （帕斯卡恒等式），代入展开。

13.1.3 概率视角#

$B_i^n(t)$ 恰好是” $n$ 次独立伯努利试验、单次成功概率为 $t$ 、恰好成功 $i$ 次”的概率。这解释了为什么 $B_i^n$ 的最大值在 $t = i/n$ 处——和二项分布的众数一致。

13.2 贝塞尔曲线的定义#

$n$ 次贝塞尔曲线由 $n+1$ 个控制点 $\mathbf{P}_0, \ldots, \mathbf{P}_n$ 定义：

\mathbf{B}(t) = \sum_{i=0}^{n} \mathbf{P}_i \, B_i^n(t), \quad t \in [0, 1]

常用次数的展开形式：

1 次（直线）： $\mathbf{B}(t) = (1-t)\mathbf{P}_0 + t\mathbf{P}_1$
2 次： $\mathbf{B}(t) = (1-t)^2\mathbf{P}_0 + 2t(1-t)\mathbf{P}_1 + t^2\mathbf{P}_2$
3 次（图形学主力）： $\mathbf{B}(t) = (1-t)^3\mathbf{P}_0 + 3t(1-t)^2\mathbf{P}_1 + 3t^2(1-t)\mathbf{P}_2 + t^3\mathbf{P}_3$

13.2.1 朴素实现（三次特化）#

GAMES101 Assignment 4 风格的直接展开：

1
cv::Point2f naive_bezier_cubic(const std::vector<cv::Point2f>& P, float t) {
2
    float u = 1.0f - t;
3
    float b0 = u * u * u;
4
    float b1 = 3.0f * u * u * t;
5
    float b2 = 3.0f * u * t * t;
6
    float b3 = t * t * t;
7
    return b0 * P[0] + b1 * P[1] + b2 * P[2] + b3 * P[3];
8
}

WARNING
朴素形式随 $n$ 增大数值稳定性下降（ $t^i$ 和 $(1-t)^{n-i}$ 数量级差很大）。实用实现请优先用下一节的 de Casteljau 算法。

13.3 de Casteljau 算法#

13.3.1 算法思想#

de Casteljau 通过递归线性插值在金字塔结构上把点 $\mathbf{B}(t)$ 算出来。它在数值上比伯恩斯坦展开稳定得多，并且几何直观。

迭代定义：

\mathbf{P}_i^{(0)} = \mathbf{P}_i, \qquad \mathbf{P}_i^{(k)} = (1-t)\mathbf{P}_i^{(k-1)} + t\mathbf{P}_{i+1}^{(k-1)}

最终结果 $\mathbf{B}(t) = \mathbf{P}_0^{(n)}$ 。

13.3.2 正确性证明#

📌 de Casteljau = 贝塞尔定义（归纳证明）#

命题： $\mathbf{P}_i^{(k)} = \sum_{j=0}^{k} \binom{k}{j} t^j (1-t)^{k-j} \mathbf{P}_{i+j}$ 。

取 $i = 0, k = n$ 即得 $\mathbf{P}_0^{(n)} = \sum_{j=0}^{n} B_j^n(t) \mathbf{P}_j = \mathbf{B}(t)$ 。

归纳基： $k = 0$ 时 $\mathbf{P}_i^{(0)} = \mathbf{P}_i$ ，平凡成立。

归纳步：设对 $k-1$ 成立。由递推式：

\mathbf{P}_i^{(k)} = (1-t)\mathbf{P}_i^{(k-1)} + t\, \mathbf{P}_{i+1}^{(k-1)}

代入归纳假设：

\mathbf{P}_i^{(k)} = (1-t)\sum_{j=0}^{k-1} \binom{k-1}{j} t^j (1-t)^{k-1-j} \mathbf{P}_{i+j} + t \sum_{j=0}^{k-1} \binom{k-1}{j} t^j (1-t)^{k-1-j} \mathbf{P}_{i+j+1}

把第二项换元 $j' = j + 1$ ，合并同类项后，利用帕斯卡恒等式 $\binom{k-1}{j} + \binom{k-1}{j-1} = \binom{k}{j}$ ：

\mathbf{P}_i^{(k)} = \sum_{j=0}^{k} \binom{k}{j} t^j (1-t)^{k-j} \mathbf{P}_{i+j} \qquad \blacksquare

13.3.3 三次贝塞尔的金字塔结构#

层级	点	表达式
0	$\mathbf{P}_0, \mathbf{P}_1, \mathbf{P}_2, \mathbf{P}_3$	原控制点
1	$\mathbf{P}_0^{(1)}, \mathbf{P}_1^{(1)}, \mathbf{P}_2^{(1)}$	相邻对的 $\text{lerp}(t)$
2	$\mathbf{P}_0^{(2)}, \mathbf{P}_1^{(2)}$	层 1 相邻对的 lerp
3	$\mathbf{P}_0^{(3)} = \mathbf{B}(t)$	最终曲线上的点

13.3.4 通用实现#

1
// 就地迭代版本，最实用，O(n²) 时间、O(n) 额外空间
2
cv::Point2f de_casteljau(std::vector<cv::Point2f> P, float t) {
3
    int n = static_cast<int>(P.size());
4
    for (int level = 1; level < n; ++level) {
5
        for (int i = 0; i + level < n; ++i) {
6
            P[i] = (1.0f - t) * P[i] + t * P[i + 1];
7
        }
8
    }
9
    return P[0];
10
}
11

12
// 教学用递归版本，清晰但有栈开销
13
cv::Point2f de_casteljau_recursive(const std::vector<cv::Point2f>& P, float t) {
14
    if (P.size() == 1) return P[0];
15
    std::vector<cv::Point2f> next(P.size() - 1);
16
    for (size_t i = 0; i + 1 < P.size(); ++i) {
17
        next[i] = (1.0f - t) * P[i] + t * P[i + 1];
18
    }
19
    return de_casteljau_recursive(next, t);
20
}

13.4 贝塞尔曲线的几何性质#

端点插值： $\mathbf{B}(0) = \mathbf{P}_0$ ， $\mathbf{B}(1) = \mathbf{P}_n$ 。

端点切向量： $\mathbf{B}'(0) = n(\mathbf{P}_1 - \mathbf{P}_0)$ ， $\mathbf{B}'(1) = n(\mathbf{P}_n - \mathbf{P}_{n-1})$ 。

这给出了设计控制点时的直观规则：要让曲线在端点有特定切线方向，只需让相邻控制点沿目标方向排布。

凸包性：曲线始终位于控制多边形的凸包内（性质 1 + 2 的直接推论）。

仿射不变性：对控制点做仿射变换 = 对曲线做同一仿射变换。这是贝塞尔曲线可以和 MVP 链自然配合的关键（先变换控制点再采样，与先采样再变换等价）。

变分递减性：任一直线与贝塞尔曲线的交点数不超过其与控制多边形的交点数——曲线”比控制多边形更平滑”的精确刻画。

13.5 贝塞尔曲线的导数#

贝塞尔曲线的导数本身是一条低一次的贝塞尔曲线，其控制点为差分：

\mathbf{B}'(t) = n \sum_{i=0}^{n-1} (\mathbf{P}_{i+1} - \mathbf{P}_i) \, B_i^{n-1}(t)

代码实现：

1
cv::Point2f bezier_derivative(const std::vector<cv::Point2f>& P, float t) {
2
    int n = static_cast<int>(P.size()) - 1;
3
    std::vector<cv::Point2f> D(n);
4
    for (int i = 0; i < n; ++i) {
5
        D[i] = static_cast<float>(n) * (P[i + 1] - P[i]);
6
    }
7
    return de_casteljau(D, t);
8
}

13.6 自适应细分#

朴素地按 $\Delta t = 0.001$ 均匀采样浪费在曲线接近直线的段上；自适应细分按曲率/平直度动态调整步长。

中点平直度判据：折半采样、检查中点到首尾连线的距离：

1
// 基于平直度递归细分
2
void adaptive_bezier(const std::vector<cv::Point2f>& P,
3
                     float t0, float t1, cv::Mat& canvas,
4
                     float tol = 0.5f, int depth = 0) {
5
    if (depth > 12) return;  // 防止病态情况下无限递归
6
    float tm = 0.5f * (t0 + t1);
7
    cv::Point2f p0 = de_casteljau(P, t0);
8
    cv::Point2f pm = de_casteljau(P, tm);
9
    cv::Point2f p1 = de_casteljau(P, t1);
10
    cv::Point2f line_mid = 0.5f * (p0 + p1);
11
    if (cv::norm(pm - line_mid) > tol) {
12
        adaptive_bezier(P, t0, tm, canvas, tol, depth + 1);
13
        adaptive_bezier(P, tm, t1, canvas, tol, depth + 1);
14
    } else {
15
        cv::line(canvas, p0, p1, cv::Scalar(0, 255, 0), 2);
16
    }
17
}

十四、B 样条与 NURBS#

14.1 为什么需要 B 样条#

贝塞尔的两个短板：

全局控制：动一个控制点，整条曲线都变——对交互建模不友好。
高次数值问题：次数随控制点数增长，既费算又不稳。

B 样条的解决方案：

分段多项式，每段只受局部几个控制点影响 → 局部支撑性
次数 $p$ 与控制点数 $n+1$ 解耦
在节点连接处自动保持 $C^{p-1}$ 连续

14.2 样条函数的定义#

设节点序列 $t_0 \leq t_1 \leq \cdots \leq t_m$ ， $p$ 次样条函数 $S(t)$ 是：

在每段 $[t_i, t_{i+1}]$ 上是次数不超过 $p$ 的多项式；
在整个定义域上具有 $C^{p-1}$ 连续性（即前 $p-1$ 阶导数在节点处连续）。

14.3 Cox-de Boor 递推公式#

14.3.1 B 样条基函数#

0 次（阶梯函数）：

N_{i,0}(t) = \begin{cases} 1, & t_i \leq t < t_{i+1} \\ 0, & \text{其他} \end{cases}

高次递推（Cox-de Boor，与 de Casteljau 形式相似）：

N_{i,p}(t) = \frac{t - t_i}{t_{i+p} - t_i} N_{i,p-1}(t) + \frac{t_{i+p+1} - t}{t_{i+p+1} - t_{i+1}} N_{i+1,p-1}(t)

WARNING
遇到重节点时分母可能为 0，约定 $0/0 := 0$ 。工程实现一定要显式处理。

14.3.2 基函数的五条性质#

非负性： $N_{i,p}(t) \geq 0$
局部支撑： $N_{i,p}(t) = 0$ 当 $t \notin [t_i, t_{i+p+1}]$ —— 这是 B 样条局部控制性的根源
单位分解： $\sum_i N_{i,p}(t) = 1$ 在基函数定义域内
节点处连续性： $C^{p-1}$ （普通节点）；重节点会降低连续性
在 $[t_i, t_{i+p+1}]$ 内恰好有一个最大值

14.4 B 样条曲线#

14.4.1 定义#

\mathbf{C}(t) = \sum_{i=0}^{n} \mathbf{P}_i \, N_{i,p}(t)

其中 $n+1$ 个控制点，节点向量长度 $m+1 = n + p + 2$ （约束： $m = n + p + 1$ ）。

14.4.2 节点向量类型#

类型	定义	特点
均匀	节点等距	基函数形状相同，平移得到
开放（clamped）	首尾节点各重复 $p+1$ 次	曲线在端点与控制多边形相切，实际建模最常用
非均匀	节点间距任意	可以精细控制局部形状

14.4.3 核心性质#

局部控制性：移动 $\mathbf{P}_i$ 只影响 $t \in [t_i, t_{i+p+1}]$ 那一段
凸包性：曲线位于相邻 $p+1$ 个控制点的局部凸包内
仿射不变性
变分递减性

14.4.4 参考实现#

1
class BSpline {
2
public:
3
    std::vector<Eigen::Vector3f> P;  // 控制点
4
    std::vector<float> U;            // 节点向量
5
    int p;                           // 次数
6

7
    // Cox-de Boor，递归实现（教学用）
8
    float N(int i, int k, float t) const {
9
        if (k == 0)
10
            return (U[i] <= t && t < U[i + 1]) ? 1.0f : 0.0f;
11
        float left = 0.0f, right = 0.0f;
12
        float d1 = U[i + k] - U[i];
13
        float d2 = U[i + k + 1] - U[i + 1];
14
        if (d1 > 1e-8f) left  = (t - U[i])         / d1 * N(i,     k - 1, t);
15
        if (d2 > 1e-8f) right = (U[i + k + 1] - t) / d2 * N(i + 1, k - 1, t);
16
        return left + right;
17
    }
18

19
    Eigen::Vector3f evaluate(float t) const {
20
        Eigen::Vector3f c = Eigen::Vector3f::Zero();
21
        for (int i = 0; i < static_cast<int>(P.size()); ++i) {
22
            c += N(i, p, t) * P[i];
23
        }
24
        return c;
25
    }
26
};

TIP
生产级实现会用 de Boor 算法（类比 de Casteljau，只在受影响的 $p+1$ 个控制点上做 $p$ 层 lerp），比直接求基函数快得多。

14.5 NURBS：有理 B 样条#

14.5.1 动机#

B 样条不能精确表示圆（以及其他二次曲线）。NURBS 通过引入控制点权重解决这一问题，CAD/工业建模里几乎全用 NURBS。

14.5.2 定义#

\mathbf{C}(t) = \frac{\sum_{i=0}^{n} w_i \mathbf{P}_i N_{i,p}(t)}{\sum_{i=0}^{n} w_i N_{i,p}(t)}

定义有理基函数：

R_{i,p}(t) = \frac{w_i N_{i,p}(t)}{\sum_{j} w_j N_{j,p}(t)}

则 $\mathbf{C}(t) = \sum_i \mathbf{P}_i R_{i,p}(t)$ ，与 B 样条形式完全一致。

14.5.3 齐次坐标视角#

把每个带权控制点升维为 $\mathbf{P}_i^w = (w_i x_i, w_i y_i, w_i z_i, w_i)$ ，则 NURBS 正是 4D 空间 B 样条的透视除法投影：

\mathbf{C}^w(t) = \sum_i \mathbf{P}_i^w N_{i,p}(t), \qquad \mathbf{C}(t) = \mathbf{C}^w(t) / w^w(t)

这个视角让 NURBS 自然与齐次坐标系统（Part 1 §1.2.3）和投影变换（Part 1 §3.1）串在一起。

14.5.4 NURBS 的三大优势#

精确表示二次曲线：圆、椭圆、抛物线、双曲线都能精确表示（有理样条的代数闭合性）
权重局部控制形状：增大 $w_i$ 把曲线向 $\mathbf{P}_i$ 拉近
投影不变性：透视投影下 NURBS 仍是 NURBS——渲染管线不会破坏它的数学结构

14.5.5 参考实现#

1
class NURBS : public BSpline {
2
public:
3
    std::vector<float> w;  // 权重，与 P 一一对应
4

5
    Eigen::Vector3f evaluate(float t) const {
6
        Eigen::Vector3f num = Eigen::Vector3f::Zero();
7
        float den = 0.0f;
8
        for (int i = 0; i < static_cast<int>(P.size()); ++i) {
9
            float b = N(i, p, t);
10
            num += w[i] * b * P[i];
11
            den += w[i] * b;
12
        }
13
        return num / den;
14
    }
15
};

14.6 贝塞尔曲面#

14.6.1 张量积定义#

双参数曲面由 $(m+1) \times (n+1)$ 的控制点网格定义：

\mathbf{S}(u, v) = \sum_{i=0}^{m} \sum_{j=0}^{n} \mathbf{P}_{ij} \, B_i^m(u) \, B_j^n(v), \quad (u, v) \in [0,1]^2

14.6.2 核心几何性质#

四角点插值： $\mathbf{S}(0,0) = \mathbf{P}_{00}$ 等四个角
边界是贝塞尔曲线：如 $\mathbf{S}(0, v) = \sum_j \mathbf{P}_{0j} B_j^n(v)$ 只由一整行/一整列控制点决定
凸包性、仿射不变性 同曲线情形

14.6.3 偏导数与法向量#

\frac{\partial \mathbf{S}}{\partial u}(u,v) = \sum_{i,j} \mathbf{P}_{ij} \frac{dB_i^m(u)}{du} B_j^n(v)

\frac{\partial \mathbf{S}}{\partial v}(u,v) = \sum_{i,j} \mathbf{P}_{ij} B_i^m(u) \frac{dB_j^n(v)}{dv}

曲面法向量由两偏导的叉积给出：

\mathbf{N}(u, v) = \frac{\partial \mathbf{S}}{\partial u} \times \frac{\partial \mathbf{S}}{\partial v}

NOTE
做光照时用的是单位法向量 $\hat{\mathbf{N}} = \mathbf{N} / \|\mathbf{N}\|$ 。如果曲面被非均匀缩放，还要按 Part 1 §2.3 的逆转置方式变换。

14.6.4 曲面细分为三角网格#

1
void tessellate_bezier_surface(const BezierSurface& s, int nu, int nv,
2
                               std::vector<Eigen::Vector3f>& out_verts,
3
                               std::vector<Eigen::Vector3f>& out_normals,
4
                               std::vector<std::array<int, 3>>& out_tris) {
5
    for (int i = 0; i <= nu; ++i) {
6
        for (int j = 0; j <= nv; ++j) {
7
            float u = static_cast<float>(i) / nu;
8
            float v = static_cast<float>(j) / nv;
9
            out_verts.push_back(s.evaluate(u, v));
10
            out_normals.push_back(s.normal(u, v));
11
        }
12
    }
13
    for (int i = 0; i < nu; ++i) {
14
        for (int j = 0; j < nv; ++j) {
15
            int i00 = i * (nv + 1) + j;
16
            int i10 = i00 + (nv + 1);
17
            int i01 = i00 + 1;
18
            int i11 = i10 + 1;
19
            out_tris.push_back({i00, i01, i10});
20
            out_tris.push_back({i01, i11, i10});
21
        }
22
    }
23
}

十五、网格几何处理#

15.1 半边数据结构（Half-Edge）#

15.1.1 为什么需要半边#

三角面 + 顶点列表（OBJ 格式那种）查询邻接关系是 $O(n)$ 的——找”和顶点 $v$ 相邻的所有面”要扫全表。半边结构用 $O(1)$ 均摊完成所有拓扑查询，是网格处理算法的标配。

15.1.2 核心约定#

每条无向边拆成两个有向半边，互为孪生（twin）
每个半边指向一个顶点（其 to），属于一条边界循环（围绕一个面）
循环方向约定：从面外看逆时针

15.1.3 数据结构#

1
struct HalfEdge {
2
    int to_vertex;   // 指向的顶点索引
3
    int face;        // 所属面（边界半边为 -1）
4
    int next;        // 同一面循环中的下一条半边
5
    int prev;        // 同一面循环中的上一条半边
6
    int twin;        // 对偶半边
7
};
8

9
struct Vertex {
10
    Eigen::Vector3f position;
11
    Eigen::Vector3f normal;
12
    int halfedge;    // 任意一条"从该顶点出发"的半边
13
};
14

15
struct Face {
16
    int halfedge;    // 该面任意一条边界半边
17
    Eigen::Vector3f normal;
18
};

15.1.4 典型 $O(1)$ 均摊查询#

1
// 遍历 v 的所有一环邻居顶点（单环遍历）
2
std::vector<int> one_ring_vertices(const HalfEdgeMesh& mesh, int v) {
3
    std::vector<int> ring;
4
    int start = mesh.vertices[v].halfedge;
5
    int he = start;
6
    do {
7
        int twin = mesh.halfedges[he].twin;
8
        ring.push_back(mesh.halfedges[twin].to_vertex);  // 注意：twin.to 才是 v 的邻居
9
        he = mesh.halfedges[twin].next;                   // 绕 v 转到下一条出边
10
    } while (he != start);
11
    return ring;
12
}

WARNING
上述循环假设 $v$ 是内部顶点。如果 $v$ 在边界上，按”孪生-next”走会掉出边界，需要特殊处理（起点选边界半边，循环退出条件改成 twin == -1）。

15.2 三角形质量度量#

网格质量直接影响后续的细分、模拟稳定性、纹理映射效果。下面三个指标最常用。

15.2.1 等周比#

Q = \frac{4\sqrt{3} \, A}{P^2}

$A$ ：三角形面积； $P$ ：周长
等边三角形时 $Q = 1$ （最优）
退化三角形 $Q \to 0$
$0 < Q \leq 1$

15.2.2 最小角度#

最小角越接近 $60°$ 越好；小于 $\sim 15°$ 的三角形叫做 sliver（针状三角形），在有限元/光照计算里会引起数值问题。

15.2.3 长宽比#

$AR = \dfrac{abc}{8(s-a)(s-b)(s-c)}$ （其中 $s$ 为半周长）。

15.2.4 实现#

1
struct TriQuality {
2
    static float isoperimetric(const Eigen::Vector3f& a,
3
                                const Eigen::Vector3f& b,
4
                                const Eigen::Vector3f& c) {
5
        float len_ab = (b - a).norm();
6
        float len_bc = (c - b).norm();
7
        float len_ca = (a - c).norm();
8
        float P = len_ab + len_bc + len_ca;
9
        if (P < 1e-8f) return 0.0f;
10
        float A = 0.5f * (b - a).cross(c - a).norm();
11
        return 4.0f * std::sqrt(3.0f) * A / (P * P);
12
    }
13

14
    static float min_angle_rad(const Eigen::Vector3f& a,
15
                                const Eigen::Vector3f& b,
16
                                const Eigen::Vector3f& c) {
17
        auto ang = [](const Eigen::Vector3f& u, const Eigen::Vector3f& v) {
18
            float cos_t = std::clamp(u.normalized().dot(v.normalized()), -1.0f, 1.0f);
19
            return std::acos(cos_t);
20
        };
21
        float A = ang(b - a, c - a);
22
        float B = ang(a - b, c - b);
23
        float C = ang(a - c, b - c);
24
        return std::min({A, B, C});
25
    }
26
};

15.3 QEM 网格简化（Quadric Error Metrics）#

Garland & Heckbert 1997 年提出的经典算法。核心思想：每次折叠最不重要的一条边、累积误差用”到邻接平面距离的平方和”度量。

15.3.1 平面方程与二次形#

平面 $ax + by + cz + d = 0$ 写成齐次形式 $\pi = (a, b, c, d)^T$ ，要求 $a^2+b^2+c^2=1$ （单位法向量）。

点 $\mathbf{p} = (x,y,z,1)^T$ 到该平面的有符号距离 = $\pi^T \mathbf{p}$ ，距离平方 = $(\pi^T \mathbf{p})^2 = \mathbf{p}^T (\pi \pi^T) \mathbf{p}$ 。

定义该平面的基本二次矩阵：

\mathbf{K}_\pi = \pi \pi^T = \begin{pmatrix} a^2 & ab & ac & ad \\ ab & b^2 & bc & bd \\ ac & bc & c^2 & cd \\ ad & bd & cd & d^2 \end{pmatrix}

15.3.2 顶点的二次误差矩阵#

顶点 $v$ 的误差矩阵是邻接面对应二次矩阵之和：

\mathbf{Q}_v = \sum_{\pi \in \text{planes}(v)} \mathbf{K}_\pi

对任意候选位置 $\mathbf{p}$ ，误差 $\Delta(\mathbf{p}) = \mathbf{p}^T \mathbf{Q}_v \mathbf{p}$ 。

15.3.3 边折叠：最优新顶点位置#

折叠边 $(v_1, v_2) \to \bar{v}$ 时：

\mathbf{Q}_{\bar{v}} = \mathbf{Q}_{v_1} + \mathbf{Q}_{v_2}

最小化 $\Delta(\mathbf{p}) = \mathbf{p}^T \mathbf{Q}_{\bar{v}} \mathbf{p}$ 关于前三个分量，令偏导为零：

\begin{pmatrix} q_{11} & q_{12} & q_{13} \\ q_{12} & q_{22} & q_{23} \\ q_{13} & q_{23} & q_{33} \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = -\begin{pmatrix} q_{14} \\ q_{24} \\ q_{34} \end{pmatrix}

若左侧矩阵奇异（退化情形），退而求其次：从 $\{v_1, v_2, (v_1+v_2)/2\}$ 三个候选中取误差最小者。

15.3.4 算法骨架#

1
// 伪代码骨架，完整实现需要堆管理 + 拓扑有效性检查
2
void qem_simplify(HalfEdgeMesh& mesh, int target_faces) {
3
    // 1. 为每个面算平面方程 -> K_pi
4
    // 2. 每个顶点累积 Q_v = ∑ K_pi
5
    // 3. 每条边算最优折叠位置 + 对应代价 Δ(p*)
6
    std::priority_queue<EdgeCost> heap;  // 最小堆
7
    for (auto& edge : mesh.edges()) {
8
        auto [pos, cost] = solve_optimal_collapse(mesh, edge);
9
        heap.push({edge.id, pos, cost});
10
    }
11
    // 4. 反复取最小代价的边折叠，直到达到目标面数
12
    while (mesh.face_count() > target_faces && !heap.empty()) {
13
        auto top = heap.top(); heap.pop();
14
        if (!mesh.is_valid_collapse(top.edge)) continue;   // 惰性删除
15
        mesh.collapse_edge(top.edge, top.position);
16
        // 更新受影响顶点的 Q_v 和邻接边的代价
17
        update_affected_edges(mesh, heap, top);
18
    }
19
}

工程上要处理的陷阱：(1) 折叠后不能产生拓扑退化（翻面、非流形边）；(2) 堆里的边项要惰性失效（用版本号或 is_valid 标志），不要真的做删除；(3) 边界边需要额外加一个”虚拟垂直墙面”保护边界形状。

15.4 Loop 细分#

15.4.1 算法概览#

针对三角网格的逼近型细分（不插值原顶点，但收敛到光滑极限曲面）。每次细分：

每条边取新中点 → 新边顶点
旧顶点按加权平均移动到新位置
每个老三角形分成 4 个新三角形

15.4.2 边中点的权重#

对内部边 $(v_1, v_2)$ ，设其两相邻三角形的另一顶点为 $v_3, v_4$ ：

\mathbf{v}_{\text{edge}} = \tfrac{3}{8}(\mathbf{v}_1 + \mathbf{v}_2) + \tfrac{1}{8}(\mathbf{v}_3 + \mathbf{v}_4)

边界边（只有一个相邻三角形）退化为：

\mathbf{v}_{\text{edge}}^{\text{bd}} = \tfrac{1}{2}(\mathbf{v}_1 + \mathbf{v}_2)

15.4.3 旧顶点的加权平均#

度数（价）为 $n$ 的内部顶点 $v$ 的新位置：

\mathbf{v}_{\text{new}} = (1 - n\beta)\mathbf{v} + \beta \sum_{i=1}^{n} \mathbf{v}_i, \qquad \beta = \begin{cases} 3/16, & n = 3 \\ 3/(8n), & n > 3 \end{cases}

边界顶点（开放网格）则用 $\mathbf{v}_{\text{new}}^{\text{bd}} = \tfrac{3}{4}\mathbf{v} + \tfrac{1}{8}(\mathbf{v}_L + \mathbf{v}_R)$ ， $\mathbf{v}_L, \mathbf{v}_R$ 为边界上的两邻居。

📌 $\beta$ 系数为何如此设计#

Loop 本人通过分析细分矩阵的特征值给出了 $\beta(n)$ ——选择让主次特征值满足 $\lambda_0 = 1$ ， $\lambda_1 = \lambda_2 = \tfrac{1}{2} + \tfrac{1}{4}\cos\tfrac{2\pi}{n}$ ，使极限曲面在普通顶点处 $C^2$ 、在奇异顶点（价不为 6 的顶点）处 $C^1$ 。细节见 Loop 1987 硕论。

15.4.4 参考实现#

1
void loop_subdivide(HalfEdgeMesh& mesh) {
2
    std::vector<Eigen::Vector3f> new_pos(mesh.vertex_count());
3
    // 1. 旧顶点新位置
4
    for (int v = 0; v < mesh.vertex_count(); ++v) {
5
        auto ring = one_ring_vertices(mesh, v);
6
        int n = static_cast<int>(ring.size());
7
        float beta = (n == 3) ? 3.0f / 16.0f : 3.0f / (8.0f * n);
8
        Eigen::Vector3f sum_neighbors = Eigen::Vector3f::Zero();
9
        for (int nb : ring) sum_neighbors += mesh.vertices[nb].position;
10
        new_pos[v] = (1.0f - n * beta) * mesh.vertices[v].position + beta * sum_neighbors;
11
    }
12
    // 2. 每条边的新顶点位置
13
    std::vector<Eigen::Vector3f> edge_pos(mesh.edge_count());
14
    for (int e = 0; e < mesh.edge_count(); ++e) {
15
        auto [v1, v2] = mesh.edge_vertices(e);
16
        auto [v3, v4] = mesh.edge_opposite_vertices(e);  // 两相邻三角形的另一点
17
        edge_pos[e] = 0.375f * (mesh.vertices[v1].position + mesh.vertices[v2].position)
18
                    + 0.125f * (mesh.vertices[v3].position + mesh.vertices[v4].position);
19
    }
20
    // 3. 重建拓扑：每个老三角形 -> 4 个新三角形
21
    mesh.rebuild_with_new_positions_and_edge_vertices(new_pos, edge_pos);
22
}

15.5 拉普拉斯平滑#

最简单的去噪 / 光顺算法。每次把顶点往邻居的重心方向拉一点：

\mathbf{v}_{\text{new}} = \mathbf{v} + \lambda \left( \frac{1}{|\mathcal{N}(v)|}\sum_{u \in \mathcal{N}(v)} \mathbf{u} - \mathbf{v} \right)

1
void laplacian_smooth(HalfEdgeMesh& mesh, float lambda = 0.5f, int iters = 10) {
2
    std::vector<Eigen::Vector3f> next(mesh.vertex_count());
3
    for (int k = 0; k < iters; ++k) {
4
        for (int v = 0; v < mesh.vertex_count(); ++v) {
5
            if (mesh.is_boundary_vertex(v)) {
6
                next[v] = mesh.vertices[v].position;  // 边界保持不动
7
                continue;
8
            }
9
            auto ring = one_ring_vertices(mesh, v);
10
            Eigen::Vector3f centroid = Eigen::Vector3f::Zero();
11
            for (int nb : ring) centroid += mesh.vertices[nb].position;
12
            centroid /= static_cast<float>(ring.size());
13
            next[v] = mesh.vertices[v].position
14
                    + lambda * (centroid - mesh.vertices[v].position);
15
        }
16
        for (int v = 0; v < mesh.vertex_count(); ++v) {
17
            mesh.vertices[v].position = next[v];
18
        }
19
    }
20
}

WARNING
朴素拉普拉斯平滑会让网格整体收缩（极限情况下塌成一点）。实用替代：Taubin λ/μ 滤波，两步交替使用正负步长，既平滑又保体积。更进阶的做法是双边保特征滤波。

小结#

主题	工具	关键性质
参数曲线	$\mathbf{C}(t)$ 、导数、曲率	图形学里曲线的通用语言；微分几何已在 Part 1 §4.1 打底
贝塞尔	伯恩斯坦基、de Casteljau	端点插值、凸包性、仿射不变性；三次最常用
B 样条	Cox-de Boor 递推、分段多项式	局部控制性、 $C^{p-1}$ 连续；次数与控制点数解耦
NURBS	有理权重 + B 样条	精确表示二次曲线；CAD/工业建模标配
半边结构	HalfEdge / Vertex / Face 三元结构	$O(1)$ 均摊拓扑查询，后续所有网格算法的基础
QEM 简化	二次误差矩阵 + 最小堆	”距离平方和”作为几何误差的代数化
Loop 细分	边中点加权 + 旧顶点加权	逼近型，收敛到光滑极限曲面
拉普拉斯平滑	邻居重心拉近	会收缩，实用时配合 Taubin λ/μ

下一部分进入 Part 4：光线追踪与全局光照——那里才是渲染方程、BRDF 完整体系、路径追踪收敛性的正式出场。

计算机图形学笔记(三)：几何建模与处理

https://kyc001.github.io/posts/计算机图形学笔记三/

作者

kyc001

发布于

2025-07-23

许可协议

CC BY-NC-SA 4.0

计算机图形学笔记(四)：光线追踪与全局光照

计算机图形学笔记(二)：光栅化渲染管线

目录#

十二、参数曲线基础#

12.1 为什么选参数形式#

12.1.1 三种曲线表示法的对比#

12.1.2 常见曲线的参数方程（速查）#

12.1.3 弧长参数化#

十三、贝塞尔曲线（canonical）#

13.1 伯恩斯坦基函数#

13.1.1 定义#

13.1.2 五条核心性质（带证明）#

13.1.3 概率视角#

13.2 贝塞尔曲线的定义#

13.2.1 朴素实现（三次特化）#

13.3 de Casteljau 算法#

13.3.1 算法思想#

13.3.2 正确性证明#

📌 de Casteljau = 贝塞尔定义（归纳证明）#

13.3.3 三次贝塞尔的金字塔结构#

13.3.4 通用实现#

13.4 贝塞尔曲线的几何性质#

13.5 贝塞尔曲线的导数#

13.6 自适应细分#

十四、B 样条与 NURBS#

14.1 为什么需要 B 样条#

14.2 样条函数的定义#

14.3 Cox-de Boor 递推公式#

14.3.1 B 样条基函数#

14.3.2 基函数的五条性质#

14.4 B 样条曲线#

14.4.1 定义#

14.4.2 节点向量类型#

14.4.3 核心性质#

14.4.4 参考实现#

14.5 NURBS：有理 B 样条#

14.5.1 动机#

14.5.2 定义#

14.5.3 齐次坐标视角#

14.5.4 NURBS 的三大优势#

14.5.5 参考实现#

14.6 贝塞尔曲面#

14.6.1 张量积定义#

14.6.2 核心几何性质#

14.6.3 偏导数与法向量#

14.6.4 曲面细分为三角网格#

十五、网格几何处理#

15.1 半边数据结构（Half-Edge）#

15.1.1 为什么需要半边#

15.1.2 核心约定#

15.1.3 数据结构#

15.1.4 典型 O(1)O(1)O(1) 均摊查询#

15.2 三角形质量度量#

15.2.1 等周比#

15.2.2 最小角度#

15.2.3 长宽比#

15.2.4 实现#

15.3 QEM 网格简化（Quadric Error Metrics）#

15.3.1 平面方程与二次形#

15.3.2 顶点的二次误差矩阵#

15.3.3 边折叠：最优新顶点位置#

15.3.4 算法骨架#

15.4 Loop 细分#

15.4.1 算法概览#

15.4.2 边中点的权重#

15.4.3 旧顶点的加权平均#

📌 β\betaβ 系数为何如此设计#

15.4.4 参考实现#

15.5 拉普拉斯平滑#

小结#

15.1.4 典型 $O(1)$ 均摊查询#

📌 $\beta$ 系数为何如此设计#