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

参数曲线的一般形式#

三维参数曲线： $\mathbf{C}(t) = \begin{pmatrix} x(t) \\ y(t) \\ z(t) \end{pmatrix}, \quad t \in [a, b]$

其中 $x(t)$、$y(t)$、$z(t)$ 是关于参数 $t$ 的连续可微函数。

参数表示相比隐式表示的优势#

1. 计算便利性：

• 参数形式：直接代入参数值即可得到曲线上的点
• 隐式形式：$F(x,y,z) = 0$ 需要求解方程

2. 动画支持： 参数 $t$ 可以直接作为时间参数，自然支持动画： $\mathbf{position}(t) = \mathbf{C}(t), \quad \mathbf{velocity}(t) = \mathbf{C}'(t)$

3. 多值函数处理： 避免了 $y = f(x)$ 形式无法表示垂直线和闭合曲线的问题。

4. 微分几何计算： 切向量、法向量、曲率等几何量的计算更加直接。

切向量与速度#

切向量（速度向量）： $\mathbf{T}(t) = \mathbf{C}'(t) = \begin{pmatrix} x'(t) \\ y'(t) \\ z'(t) \end{pmatrix}$

单位切向量： $\hat{\mathbf{T}}(t) = \frac{\mathbf{T}(t)}{\|\mathbf{T}(t)\|} = \frac{\mathbf{C}'(t)}{\|\mathbf{C}'(t)\|}$

曲率的数学定义#

曲率的几何定义： 曲率 $\kappa(t)$ 描述曲线在某点处偏离直线的程度： $\kappa(t) = \frac{\|\mathbf{T}'(t)\|}{\|\mathbf{C}'(t)\|} = \frac{\|\mathbf{C}'(t) \times \mathbf{C}''(t)\|}{\|\mathbf{C}'(t)\|^3}$

平面曲线的曲率公式： 对于平面曲线 $\mathbf{C}(t) = (x(t), y(t))$： $\kappa(t) = \frac{|x'(t)y''(t) - y'(t)x''(t)|}{(x'(t)^2 + y'(t)^2)^{3/2}}$

弧长参数化#

弧长函数： $s(t) = \int_a^t \|\mathbf{C}'(\tau)\| d\tau$

弧长参数化的优势： 当曲线以弧长为参数时，$\|\mathbf{C}'(s)\| = 1$，简化了几何计算： $\kappa(s) = \|\mathbf{C}''(s)\|$

代码实现：

1
float calculate_curvature_2d(float t,
2
                            std::function<Vector2f(float)> curve,
3
                            std::function<Vector2f(float)> first_derivative,
4
                            std::function<Vector2f(float)> second_derivative) {
5
    Vector2f first = first_derivative(t);
6
    Vector2f second = second_derivative(t);
7

8
    float numerator = std::abs(first.x() * second.y() - first.y() * second.x());
9
    float denominator = std::pow(first.squaredNorm(), 1.5f);
10

11
    return numerator / denominator;
12
}

线性插值的数学形式#

两点间直线的参数方程： $\mathbf{L}(t) = \mathbf{P}_0 + t(\mathbf{P}_1 - \mathbf{P}_0) = (1-t)\mathbf{P}_0 + t\mathbf{P}_1, \quad t \in [0, 1]$

几何意义：

• $t = 0$：位于点 $\mathbf{P}_0$
• $t = 1$：位于点 $\mathbf{P}_1$
• $0 < t < 1$：位于线段 $\mathbf{P}_0\mathbf{P}_1$ 内部

微分性质：

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

工程实现#

1
class ParametricLine {
2
private:
3
    Vector3f p0, p1;
4

5
public:
6
    ParametricLine(const Vector3f& start, const Vector3f& end) : p0(start), p1(end) {}
7

8
    Vector3f evaluate(float t) const {
9
        return (1.0f - t) * p0 + t * p1;
10
    }
11

12
    Vector3f tangent() const {
13
        return (p1 - p0).normalized();
14
    }
15

16
    float length() const {
17
        return (p1 - p0).norm();
18
    }
19
};

圆的参数方程#

标准圆的参数形式： $\mathbf{C}(t) = \mathbf{center} + r\begin{pmatrix} \cos(t) \\ \sin(t) \end{pmatrix}, \quad t \in [0, 2\pi]$

微分性质：

• 切向量：$\mathbf{C}'(t) = r\begin{pmatrix} -\sin(t) \\ \cos(t) \end{pmatrix}$
• 速度大小：$\|\mathbf{C}'(t)\| = r$
• 曲率：$\kappa = \frac{1}{r}$（常曲率）

椭圆的参数方程#

标准椭圆： $\mathbf{E}(t) = \mathbf{center} + \begin{pmatrix} a\cos(t) \\ b\sin(t) \end{pmatrix}, \quad t \in [0, 2\pi]$

其中 $a$ 和 $b$ 分别是长轴和短轴的半长度。

椭圆的曲率： $\kappa(t) = \frac{ab}{(a^2\sin^2(t) + b^2\cos^2(t))^{3/2}}$

特殊点的曲率：

• 长轴端点：$\kappa = \frac{b^2}{a^3}$
• 短轴端点：$\kappa = \frac{a^2}{b^3}$

三维空间中的圆#

任意平面上的圆： $\mathbf{C}(t) = \mathbf{center} + r(\cos(t)\mathbf{u} + \sin(t)\mathbf{v})$

其中 $\mathbf{u}$ 和 $\mathbf{v}$ 是圆所在平面的两个正交单位向量。

伯恩斯坦多项式的数学定义#

严格定义： $n$ 次伯恩斯坦基函数定义为： $B_i^n(t) = \binom{n}{i} t^i (1-t)^{n-i}, \quad i = 0, 1, \ldots, n$

其中二项式系数： $\binom{n}{i} = \frac{n!}{i!(n-i)!}$

基本性质的数学证明#

1. 非负性：$B_i^n(t) \geq 0$ 对所有 $t \in [0,1]$

• 证明：$\binom{n}{i} \geq 0$，$t^i \geq 0$，$(1-t)^{n-i} \geq 0$

2. 归一性（权重和为1）：$\sum_{i=0}^{n} B_i^n(t) = 1$

• 证明：利用二项式定理 $\sum_{i=0}^{n} B_i^n(t) = \sum_{i=0}^{n} \binom{n}{i} t^i (1-t)^{n-i} = (t + (1-t))^n = 1$

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}$（Kronecker delta）
• $B_i^n(1) = \delta_{in}$

5. 递推关系： $B_i^n(t) = (1-t)B_i^{n-1}(t) + tB_{i-1}^{n-1}(t)$

这个递推关系是de Casteljau算法的数学基础。

伯恩斯坦多项式的几何意义#

概率解释： $B_i^n(t)$ 可以解释为在 $n$ 次独立伯努利试验中，成功概率为 $t$ 时，恰好成功 $i$ 次的概率。

权重函数性质：

• 在 $t = 0$ 时，只有 $B_0^n(0) = 1$，其他为0
• 在 $t = 1$ 时，只有 $B_n^n(1) = 1$，其他为0
• 在 $t = i/n$ 时，$B_i^n(t)$ 达到最大值

de Casteljau算法的几何原理#

核心思想： de Casteljau算法通过递归的线性插值来计算贝塞尔曲线上的点，这种方法在数值上比直接计算伯恩斯坦多项式更稳定。

数学表述： 设控制点为 $\mathbf{P}_0, \mathbf{P}_1, \ldots, \mathbf{P}_n$，定义递归序列： $\mathbf{P}_i^{(0)} = \mathbf{P}_i, \quad i = 0, 1, \ldots, n$ $\mathbf{P}_i^{(k)} = (1-t)\mathbf{P}_i^{(k-1)} + t\mathbf{P}_{i+1}^{(k-1)}, \quad k = 1, 2, \ldots, n; \quad i = 0, 1, \ldots, n-k$

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

算法的数学证明#

定理：de Casteljau算法计算的结果等于贝塞尔曲线的定义式。

证明思路： 利用数学归纳法证明： $\mathbf{P}_i^{(k)} = \sum_{j=0}^{k} \binom{k}{j} t^j (1-t)^{k-j} \mathbf{P}_{i+j}$

基础步骤：$k=0$ 时显然成立。

归纳步骤：假设对 $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+1+j}$

化简得到： $\mathbf{P}_i^{(k)} = \sum_{j=0}^{k} \binom{k}{j} t^j (1-t)^{k-j} \mathbf{P}_{i+j}$

算法的几何解释#

三次贝塞尔曲线示例： 对于控制点 $\mathbf{P}_0, \mathbf{P}_1, \mathbf{P}_2, \mathbf{P}_3$：

第一层插值：

• $\mathbf{P}_0^{(1)} = (1-t)\mathbf{P}_0 + t\mathbf{P}_1$
• $\mathbf{P}_1^{(1)} = (1-t)\mathbf{P}_1 + t\mathbf{P}_2$
• $\mathbf{P}_2^{(1)} = (1-t)\mathbf{P}_2 + t\mathbf{P}_3$

第二层插值：

• $\mathbf{P}_0^{(2)} = (1-t)\mathbf{P}_0^{(1)} + t\mathbf{P}_1^{(1)}$
• $\mathbf{P}_1^{(2)} = (1-t)\mathbf{P}_1^{(1)} + t\mathbf{P}_2^{(1)}$

第三层插值：

• $\mathbf{P}_0^{(3)} = (1-t)\mathbf{P}_0^{(2)} + t\mathbf{P}_1^{(2)}$

最终 $\mathbf{B}(t) = \mathbf{P}_0^{(3)}$

样条函数的严格定义#

样条函数： 设节点序列 $t_0 < t_1 < \cdots < t_m$，$p$ 次样条函数 $S(t)$ 是满足以下条件的分段多项式：

1. 在每个区间 $[t_i, t_{i+1}]$ 上，$S(t)$ 是次数不超过 $p$ 的多项式
2. $S(t)$ 在整个定义域上具有 $C^{p-1}$ 连续性

B样条基函数的递推定义#

0次B样条基函数：

基函数定义：

当 $t_i \leq t < t_{i+1}$ 时：$N_{i,0}(t) = 1$

其他情况：$N_{i,0}(t) = 0$

高次B样条基函数（Cox-de Boor递推公式）： $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)$

约定：当分母为零时，对应项为零，即 $\frac{0}{0} = 0$。

B样条基函数的重要性质#

1. 非负性：$N_{i,p}(t) \geq 0$ 对所有 $t$

2. 局部支撑性：$N_{i,p}(t) = 0$ 当 $t \notin [t_i, t_{i+p+1}]$

3. 权重和为1：$\sum_{i} N_{i,p}(t) = 1$ 对所有 $t$

4. 连续性：$N_{i,p}(t)$ 在节点处具有 $C^{p-1}$ 连续性

B样条曲线的定义#

数学表达式： $\mathbf{C}(t) = \sum_{i=0}^{n} \mathbf{P}_i N_{i,p}(t)$

其中：

• $\mathbf{P}_i$ 是控制点
• $N_{i,p}(t)$ 是 $p$ 次B样条基函数
• $n+1$ 是控制点的数量

节点向量的作用#

节点向量：$\mathbf{T} = \{t_0, t_1, \ldots, t_{m}\}$，其中 $m = n + p + 1$

节点向量的分类：

1. 均匀节点向量：节点等间距分布
2. 非均匀节点向量：节点间距不等
3. 开放节点向量：首末节点重复度为 $p+1$

B样条曲线的性质#

1. 凸包性质：曲线位于控制点的凸包内

2. 局部控制性：移动控制点 $\mathbf{P}_i$ 只影响参数区间 $[t_i, t_{i+p+1}]$ 内的曲线

3. 变分递减性：曲线与任意直线的交点数不超过控制多边形与该直线的交点数

4. 仿射不变性：仿射变换可直接作用于控制点

代码实现：

1
class BSpline {
2
private:
3
    std::vector<Vector3f> control_points;
4
    std::vector<float> knot_vector;
5
    int degree;
6

7
public:
8
    float basis_function(int i, int p, float t) {
9
        if (p == 0) {
10
            return (t >= knot_vector[i] && t < knot_vector[i + 1]) ? 1.0f : 0.0f;
11
        }
12

13
        float left_coeff = 0.0f, right_coeff = 0.0f;
14

15
        if (knot_vector[i + p] != knot_vector[i]) {
16
            left_coeff = (t - knot_vector[i]) / (knot_vector[i + p] - knot_vector[i]);
17
        }
18

19
        if (knot_vector[i + p + 1] != knot_vector[i + 1]) {
20
            right_coeff = (knot_vector[i + p + 1] - t) / (knot_vector[i + p + 1] - knot_vector[i + 1]);
21
        }
22

23
        return left_coeff * basis_function(i, p - 1, t) +
24
               right_coeff * basis_function(i + 1, p - 1, t);
25
    }
26

27
    Vector3f evaluate(float t) {
28
        Vector3f point = Vector3f::Zero();
29

30
        for (int i = 0; i < control_points.size(); ++i) {
31
            point += control_points[i] * basis_function(i, degree, t);
32
        }
33

34
        return point;
35
    }
36
};

NURBS的定义#

非均匀有理B样条（NURBS）： $\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)}$

其中 $w_i$ 是控制点 $\mathbf{P}_i$ 的权重。

有理基函数#

有理基函数定义： $R_{i,p}(t) = \frac{w_i N_{i,p}(t)}{\sum_{j=0}^{n} w_j N_{j,p}(t)}$

NURBS曲线的简化表示： $\mathbf{C}(t) = \sum_{i=0}^{n} \mathbf{P}_i R_{i,p}(t)$

NURBS的几何意义#

齐次坐标表示： NURBS可以看作是四维空间中B样条曲线的透视投影： $\mathbf{C}^w(t) = \sum_{i=0}^{n} \mathbf{P}_i^w N_{i,p}(t)$

其中 $\mathbf{P}_i^w = (w_i x_i, w_i y_i, w_i z_i, w_i)$ 是齐次控制点。

透视投影： $\mathbf{C}(t) = \frac{(\mathbf{C}^w_x(t), \mathbf{C}^w_y(t), \mathbf{C}^w_z(t))}{\mathbf{C}^w_w(t)}$

NURBS的优势#

1. 精确表示圆锥曲线：通过适当选择权重，可以精确表示圆、椭圆、抛物线、双曲线

2. 局部权重控制：改变权重 $w_i$ 可以局部调整曲线形状

3. 投影不变性：透视投影下NURBS仍为NURBS

代码实现：

1
class NURBS : public BSpline {
2
private:
3
    std::vector<float> weights;
4

5
public:
6
    Vector3f evaluate(float t) override {
7
        Vector3f numerator = Vector3f::Zero();
8
        float denominator = 0.0f;
9

10
        for (int i = 0; i < control_points.size(); ++i) {
11
            float basis = basis_function(i, degree, t);
12
            numerator += weights[i] * control_points[i] * basis;
13
            denominator += weights[i] * basis;
14
        }
15

16
        return numerator / denominator;
17
    }
18
};

贝塞尔曲面的定义#

双参数贝塞尔曲面： $\mathbf{S}(u,v) = \sum_{i=0}^{m} \sum_{j=0}^{n} \mathbf{P}_{ij} B_i^m(u) B_j^n(v)$

其中：

• $\mathbf{P}_{ij}$ 是 $(m+1) \times (n+1)$ 控制点网格
• $B_i^m(u)$ 和 $B_j^n(v)$ 是伯恩斯坦基函数
• $(u,v) \in [0,1] \times [0,1]$ 是参数域

曲面的几何性质#

边界曲线： 贝塞尔曲面的四条边界都是贝塞尔曲线：

四条边界曲线：

左边界（$u=0$）： $\mathbf{S}(0,v) = \sum_{j=0}^{n} \mathbf{P}_{0j} B_j^n(v)$

右边界（$u=1$）： $\mathbf{S}(1,v) = \sum_{j=0}^{n} \mathbf{P}_{mj} B_j^n(v)$

下边界（$v=0$）： $\mathbf{S}(u,0) = \sum_{i=0}^{m} \mathbf{P}_{i0} B_i^m(u)$

上边界（$v=1$）： $\mathbf{S}(u,1) = \sum_{i=0}^{m} \mathbf{P}_{in} B_i^m(u)$

角点插值： $\mathbf{S}(0,0) = \mathbf{P}_{00}, \quad \mathbf{S}(1,0) = \mathbf{P}_{m0}, \quad \mathbf{S}(0,1) = \mathbf{P}_{0n}, \quad \mathbf{S}(1,1) = \mathbf{P}_{mn}$

曲面法向量计算#

偏导数： $\frac{\partial \mathbf{S}}{\partial u} = \sum_{i=0}^{m} \sum_{j=0}^{n} \mathbf{P}_{ij} \frac{dB_i^m(u)}{du} B_j^n(v)$

$\frac{\partial \mathbf{S}}{\partial v} = \sum_{i=0}^{m} \sum_{j=0}^{n} \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}$

单位法向量： $\hat{\mathbf{N}}(u,v) = \frac{\mathbf{N}(u,v)}{\|\mathbf{N}(u,v)\|}$

代码实现：

1
class BezierSurface {
2
private:
3
    std::vector<std::vector<Vector3f>> control_points;  // m+1 x n+1 网格
4
    int m, n;  // 度数
5

6
public:
7
    Vector3f evaluate(float u, float v) {
8
        Vector3f point = Vector3f::Zero();
9

10
        for (int i = 0; i <= m; ++i) {
11
            for (int j = 0; j <= n; ++j) {
12
                float basis_u = bernstein_polynomial(i, m, u);
13
                float basis_v = bernstein_polynomial(j, n, v);
14
                point += control_points[i][j] * basis_u * basis_v;
15
            }
16
        }
17

18
        return point;
19
    }
20

21
    Vector3f normal(float u, float v) {
22
        Vector3f du = partial_derivative_u(u, v);
23
        Vector3f dv = partial_derivative_v(u, v);
24
        return du.cross(dv).normalized();
25
    }
26

27
private:
28
    float bernstein_polynomial(int i, int n, float t) {
29
        return binomial_coefficient(n, i) *
30
               std::pow(t, i) * std::pow(1 - t, n - i);
31
    }
32
};

三角形质量度量#

等周比（Isoperimetric Ratio）： $Q = \frac{4\sqrt{3} \cdot A}{P^2}$

其中：

• $A$ 是三角形面积
• $P$ 是三角形周长

数学性质：

• $Q = 1$ 对于等边三角形（最优）
• $Q \to 0$ 对于退化三角形
• $0 < Q \leq 1$ 对于所有有效三角形

面积计算： 对于顶点 $\mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3$： $A = \frac{1}{2}\|(\mathbf{v}_2 - \mathbf{v}_1) \times (\mathbf{v}_3 - \mathbf{v}_1)\|$

周长计算： $P = \|\mathbf{v}_2 - \mathbf{v}_1\| + \|\mathbf{v}_3 - \mathbf{v}_2\| + \|\mathbf{v}_1 - \mathbf{v}_3\|$

其他质量度量#

长宽比（Aspect Ratio）： $AR = \frac{a \cdot b \cdot c}{8(s-a)(s-b)(s-c)}$

其中 $s = \frac{a+b+c}{2}$ 是半周长。

最小角度： $\theta_{min} = \min\{\arccos\frac{b^2+c^2-a^2}{2bc}, \arccos\frac{a^2+c^2-b^2}{2ac}, \arccos\frac{a^2+b^2-c^2}{2ab}\}$

工程实现#

1
class TriangleQuality {
2
public:
3
    static float isoperimetric_ratio(const Vector3f& v1, const Vector3f& v2, const Vector3f& v3) {
4
        Vector3f e1 = v2 - v1;
5
        Vector3f e2 = v3 - v2;
6
        Vector3f e3 = v1 - v3;
7

8
        float a = e1.norm();
9
        float b = e2.norm();
10
        float c = e3.norm();
11

12
        float area = 0.5f * e1.cross(-e3).norm();  // 使用叉积计算面积
13
        float perimeter = a + b + c;
14

15
        if (perimeter < 1e-8f) return 0.0f;  // 避免除零
16

17
        return 4.0f * std::sqrt(3.0f) * area / (perimeter * perimeter);
18
    }
19

20
    static float min_angle(const Vector3f& v1, const Vector3f& v2, const Vector3f& v3) {
21
        Vector3f e1 = (v2 - v1).normalized();
22
        Vector3f e2 = (v3 - v2).normalized();
23
        Vector3f e3 = (v1 - v3).normalized();
24

25
        float angle1 = std::acos(std::clamp(e1.dot(-e3), -1.0f, 1.0f));
26
        float angle2 = std::acos(std::clamp((-e1).dot(e2), -1.0f, 1.0f));
27
        float angle3 = std::acos(std::clamp((-e2).dot(e3), -1.0f, 1.0f));
28

29
        return std::min({angle1, angle2, angle3});
30
    }
31
};

网格统计信息：