第21章：隐式建模

21.1 隐式函数与骨架图元

21.1.1 什么是隐式曲面

前述章节的几何表示主要是显式的：参数曲面 p(u,v) 或三角形网格，直接给出曲面上的每个点。而隐式表示（Implicit Representation）截然不同——它用一个标量函数 f : ℝ³ → ℝ 的零等值面来定义曲面。

══════════════════════════════════════════════
  隐式曲面 (Implicit Surface) — 正式定义
══════════════════════════════════════════════

// 标量场 f : ℝ³ → ℝ 在全空间每点赋一个实数值
// 零等值面（Zero Isosurface）即隐式曲面：

S = { p ∈ ℝ³ | f(p) = 0 }

// 空间三区域划分：
f(p) < 0  →  p 位于物体内部 (Interior)
f(p) = 0  →  p 位于物体曲面上 (Surface)
f(p) > 0  →  p 位于物体外部 (Exterior)

// 梯度给出曲面法向量方向：
∇f(p) = ( ∂f/∂x , ∂f/∂y , ∂f/∂z )ᵀ

// 单位法向量 (假定 ||∇f|| ≠ 0)：
n̂(p) = ∇f(p) / ||∇f(p)||

// 隐式曲面的数学优势：
//   • 内/外判定 O(1)：仅需计算 f(p) 的符号
//   • CSG 布尔操作 O(1)：min/max 即得结果
//   • 拓扑变化自然：场值连续变化 → 曲面自动融合/分裂
//   • 碰撞检测极速：有 SDF 时 f(p) 即到曲面距离

曲面将空间划分为三个区域：
f(p) < 0 — 内部（物体内部）
f(p) = 0 — 曲面（零等值面）
f(p) > 0 — 外部（物体外部）

这一约定虽非强制，但在 CSG 操作中极为便利：并取 min 等价于"内部区域的并集"，交取 max 等价于"内部区域的交集"。

梯度（Gradient）给出了曲面在每个点的朝向：

向量 ∇f(p) 在 f(p)=0 上垂直于曲面，提供了该点的法向量（未归一化时）。与参数表示不同，隐式表示不需要显式存储法线——它在任何点都可通过偏导数计算得到。

21.1.2 隐式表示 vs 显式表示：精度对比

21.1.3 骨架图元：Blob、椭球与圆柱

构造隐式曲面时，在空间中放置若干骨架图元（Skeletal Primitives），每个贡献一个标量场。图元的场函数在骨架中心处取最大值，随距离单调衰减至零。

点骨架（Blob）—— Blinn 高斯模型

最基础的图元是点骨架（Blob），由位置 c 定义。1982年，Blinn 受电子密度云启发，提出用高斯函数定义场：

引入阈值 T 后定义等值面：

令 F(p) = 0 解出等值面的有效半径：

Blinn 模型的梯度：

与参数 b 成正比——b 越大衰减越快，梯度越陡峭。这对数值稳定性有直接影响：b 过大时等值面附近场值急剧变化，使 Marching Cubes 的线性插值假设失效，需极小体素才能准确捕获。

Wyvill 六次多项式

Wyvill 兄弟（1986）提出用紧致支撑（Compact Support）多项式替代高斯函数，彻底解决效率问题：

Wyvill 多项式的展开形式（令 R=1）——多项式系数的精确值决定函数在支撑边界的各阶光滑性：

// Wyvill 六次多项式 — 紧致支撑 C² 场函数
// 原始形式：f(r) = (1 - r²)³     (0 ≤ r ≤ 1)
// 展开形式揭示各阶导数连续性：

f(r) = 1 - 3r² + 3r⁴ - r⁶

// 另一种常用归一化形式 (Wyvill et al. 1986)：
f(r) = 1 - (22/9)r² + (17/9)r⁴ - (4/9)r⁶

// 等价变形：f(r) = 4r⁶/9 - 17r⁴/9 + 22r²/9  (归一化形式，值域 [0,1])

// 各阶导数在 r=0 处的值：
f(0) = 1          // 峰值
f'(0) = 0          // 平滑顶部 (一阶导为零)
f''(0) = -6/R²     // 有限曲率

// 各阶导数在 r=R 处的值 (紧致支撑的关键)：
f(R) = 0          // 边界值为零
f'(R) = 0          // 一阶导连续为零 → C¹
f''(R) = 0         // 二阶导连续为零 → C² (视觉无折痕)

关键性质——在边界 r=R 处的连续性分析：

函数在中心 r=0 处：f(0)=1（峰值），f'(0)=0（平滑顶部），f''(0)=−6/R²。

椭球体与圆柱体

椭球体（Ellipsoid）通过变换矩阵将球拉伸为椭球。设变换 M = S · R（缩放 S × 旋转 R），局部坐标距离：

若 M = diag(1, 1, 1)，退化为球体；若 M = diag(a, b, c) 则为轴对齐椭球，半轴长度分别为 a、b、c。

圆柱体（Cylinder）距离为点到中心线段的垂直距离。设线段端点 a, b：

沿轴线方向的场为常值（无限长圆柱），如需截断可在两端加半球帽（capsule 形状）。

骨架图元场函数速查表

21.1.4 组合多个图元：求和混合

多个骨架图元的场可通过求和混合组合为一个总场：

这是隐式建模中最简单也最经典的组合方式。当两个 blob 相互靠近时，它们的场在空间叠加：重叠区域的场值升高（因为两正值相加），导致零等值面向外"鼓胀"，形成平滑的哑铃形状（Surface Fusion）。

Blinn 原始公式（1982）直接使用了这种加权求和形式：

求和混合的直观理解：每个图元在其支撑域内贡献"物质密度"，空间中某点的总密度超过阈值 T 即视为物体内部。T 越大物体越小（等值面收缩），T 越小物体越大（等值面膨胀）。

21.2 距离场与等值面

21.2.1 有符号距离函数（SDF）的精确定义

有符号距离函数（Signed Distance Function，SDF）是隐式建模中一种特殊的标量场——它不仅给出到曲面的距离，还保留内/外符号信息。

形式上，定义曲面 S 将空间分为内部 Ω− 和外部 Ω₄。最近点映射 π(p) = argmin˰˲˯˭ˬ˭ ||p−q||，则：

SDF 的核心性质：Eikonal 方程

距离函数的梯度为单位向量（除不可微的脊线外）：

══════════════════════════════════════════════
  有符号距离函数 (SDF) — Eikonal 方程与性质
══════════════════════════════════════════════

// SDF 定义：标量场 f(p) 给出点 p 到曲面 ∂Ω 的最短距离
//   符号为正表示 p 在物体外部，符号为负表示在内部

f(p) = sign(p) · min_{q ∈ ∂Ω} ||p − q||

// Eikonal 方程 — SDF 的核心性质 (几乎处处成立)：
||∇f(p)|| = 1

// 推论 1：光线步进安全步长 → step = |f(p)| 不越过曲面
// 推论 2：法线无需归一化 → n̂(p) = ∇f(p) (已为单位向量)
// 推论 3：厚度偏移 → f(p)−d = 0 是原曲面沿法线外移 d 的等距面

// ─── CSG 操作 (保持 SDF 性质，除操作交界处外) ───
// 并 (Union)：两点取更近者 (更靠内部的 min)
f_∪ = min( f₁, f₂ )

// 交 (Intersection)：两点取更远者 (更靠外部的 max)
f_∩ = max( f₁, f₂ )

// 差 (Subtraction)：内部取反 → A 中挖去 B
f_− = max( f₁, −f₂ )

// 注：SDF 在曲面"脊线"处不可微 (距离函数不光滑)
//    脊线是空间中到曲面有 ≥2 个最近点的集合 (Medial Axis)

证明概要：在光滑区域，沿梯度方向移动 ε 使函数值变化 ε||∇f||。但沿法线方向移动 ε 恰好使距离变化 ε（因为这是到曲面的直线方向）。因此 ||∇f|| = 1。

这一性质的三个重要推论：
(1) 光线步进安全步长：沿光线方向，f(p) 的绝对值即为到曲面的最小距离——可以直接安全步进 ||f(p)|| 而不越过曲面。
(2) 法线无需归一化：∇f(p) 本身就是单位法向量（当然仅在微分存在处）。
(3) 厚度偏移：f(p) − d = 0 定义为原曲面沿法线偏移 d 的等距面。

21.2.2 基本 SDF 公式集

══════════════════════════════════════════════
  基本 SDF 图元 — 完整公式集
══════════════════════════════════════════════

// ─── 球体 (Sphere) ───
f_sphere(p, c, R) = ||p − c|| − R
// ∇f = (p−c)/||p−c|| → ||∇f|| = 1 ✓

// ─── 长方体 (Box, 半边长 b=(b_x,b_y,b_z)) ───
q = |p − c| − b
f_box(p) = ||max(q, 0)|| + min(max(q_x, q_y, q_z), 0)
// 第一项：外部点到最近边的距离；第二项：内部负距离

// ─── 圆环面 (Torus, 大径 R, 小径 r) ───
f_torus(p) = || ( ||(p_x, p_z)|| − R ,  p_y ) || − r
// 先将 p 投影到 xz 平面求到圆心距离，再在 2D (ρ,y) 中求到圆环距离

// ─── 胶囊体 (Capsule, 端点 a,b, 半径 R) ───
t = clamp( (p−a)·(b−a) / ||b−a||² , 0, 1 )
f_capsule(p) = || p − (a + t·(b−a)) || − R

// ─── 圆柱体 (无限长沿 y 轴, 半径 R) ───
f_cylinder(p) = √(p_x² + p_z²) − R

21.2.3 SDF 上的 CSG 操作

SDF 支持精确的构造实体几何（CSG）操作，通过简单的最值函数即可实现：

21.2.4 等值面与水平集

等值面（Isosurface）是标量场 f 中所有满足某个常数值的点的集合。零等值面 f(p)=0 是最常用的曲面定义，但任意等值面 f(p)=c 同样有意义。

改变水平值产生"嵌套"曲面族——在 SDF 中，f(p)=δ 是原曲面沿法线方向外偏移 δ 的等距面。这使得隐式表示特别适合厚度壳体建模：取两等值面 f(p)=+δ 和 f(p)=−δ 之间的区域即为厚度 2δ 的壳体。

法线计算：曲面上任一点 p（满足 f(p)=0）的单位法向量由梯度给出：

══════════════════════════════════════════════
  SDF 几何变换与梯度计算
══════════════════════════════════════════════

// ─── 刚体变换 (保持距离场性质) ───
// 平移：对输入点施加逆平移
f_trans(p) = f(p − t)              t 是平移向量

// 旋转：对输入点施加逆旋转 (R⁻¹ = Rᵀ 对正交矩阵)
f_rot(p) = f( R⁻¹ · p )            R 是 3×3 旋转矩阵

// 均匀缩放 (保持 ||∇f||=1 性质!)：
f_scale(p) = s · f(p/s)            s>0 是均匀缩放因子

// 验证：|∇f_scale| = s·|∇f(p/s)|·(1/s) = |∇f(p/s)| = 1 ✓

// ─── 梯度计算 (有限差分) ───
// 解析法：对 Wyvill f(r) 链式求导
∇f(p) = f'(r) · (∇r)
f'(r) = −6r/R² · (1−r²/R²)²       Wyvill 一阶导解析式

// 数值法：中心差分 (用于复杂混合场的梯度)
∂f/∂x ≈ ( f(p_x+h, p_y, p_z) − f(p_x−h, p_y, p_z) ) / (2h)
// h = max(ε, ε·|p|) 自适应步长防止浮点下溢 (ε≈10⁻⁶)

在 Blob-tree 混合场中，梯度可通过解析求导或有限差分 (±h 采样) 近似。解析求导更精确且效率更高，是现代实现的标准做法。

21.3 空间划分

21.3.1 为什么需要空间划分

考虑一个含 N=100,000 个 Wyvill blob 的场景。对空间任一点 p，直接计算：

Marching Cubes 需对网格上 M 个体素（通常 M=10&sup6;~10⁷）各求值 8 个顶点。直接暴力求值总代价 O(N·M) ≈ 10⁶−10¹¹，完全不可行。

空间划分的核心思想：只访问对 p 有非零贡献的图元。利用 Wyvill 函数的紧致支撑特性，将图元按空间位置组织，查询时仅检索覆盖 p 的局部图元。

21.3.2 均匀网格（Uniform Grid）

将包围盒 [xᵢᵗᵗ, xᵢᵜᵤ] × [yᵢᵗᵗ, yᵢᵜᵤ] × [zᵢᵗᵗ, zᵢᵜᵤ] 均匀剖分为 n₆ × n₇ × nₙ 个立方体单元，边长 Δ。每个图元按其支撑域（以骨架中心 c 为球心、半径 R 的球）覆盖的单元集合注册到各单元。

查询点 p 时，O(1) 定位到所在单元，仅求值该单元内存储的图元。设每个单元平均存储 K 个图元，则单点求值 O(K)。

优缺点：实现简单，查询快速。但若非均匀场景（局部 blob 密集、其他区域稀疏），网格边长难以统一选定：太大则密集单元 K 过大，太小则整体内存 O(n³) 爆炸。

21.3.3 八叉树（Octree）自适应空间划分

八叉树（Octree）自适应调整分辨率，在稀疏区用粗粒度、在密集区用细粒度。

构建算法

分裂准则

决定是否继续细分有两种主要策略：

查询点 p 时从根节点递归下降：在每层判断 p 位于哪个子节点，直到叶子节点，累计所有经过节点的图元（或仅求值叶子层图元，取决于实现）。查询复杂度 O(log N · K)，其中 N 为节点数。

性能对比

21.3.4 自适应求值策略

自适应求值（Adaptive Evaluation）在零等值面附近（|f(p)| < ε）用高分辨率，远处用粗粒度近似。这与自适应 Marching Cubes（见 21.5 节）在概念上一致，但作用于求值而非多边形化阶段。

实现要点：
(1) 深度上限：限制八叉树最大深度（通常 8–12 层），防止单个体素过小导致浮点精度问题。
(2) 惰性细化：初始构建至中等深度，查询时若某节点邻接零等值面则按需细化。
(3) 缓存友好：将节点数据按 Morton 码（Z-order curve）排序存储，使空间上靠近的节点在内存中也靠近，提升缓存命中率。

══════════════════════════════════════════════
  Morton 码 (Z-order) — 八叉树缓存优化
══════════════════════════════════════════════

// Morton 码：将 3D 坐标 (x,y,z) 的二进制位交错排列
//   使空间邻近的点在内存中也邻近 → 缓存命中率 ↑

// 编码 (interleave bits of x,y,z)：
morton(x,y,z) = Σ_{i=0}^{31} (bit_i(x)·4 + bit_i(y)·2 + bit_i(z)) · 8^i
// 示例：x=1(001), y=3(011), z=2(010) → Morton = 0_1_0_0_1_1_0_0_1 = 73

// 高效实现 (位运算扩展每个坐标到三倍宽度再交错)：
expand_bits(v):  // 将 10 位整数扩展为 30 位 (每 bit 后插 2 个零)
  v = (v | (v << 16)) & 0x030000FF
  v = (v | (v <<  8)) & 0x0300F00F
  v = (v | (v <<  4)) & 0x030C30C3
  v = (v | (v <<  2)) & 0x09249249
  return v

morton_encode(x,y,z):
  return expand_bits(x) | (expand_bits(y) << 1) | (expand_bits(z) << 2)

21.4 更多混合方式

══════════════════════════════════════════════
  混合算子完整公式汇编 (Blending Operators)
══════════════════════════════════════════════

// ─── ① 求和混合 (Sum Blend) ───
F_sum = (Σᵢ fᵢ(p)) − T
// 最简单但产生"涨大效应" (Bulging)

// ─── ② Ricci 混合 (幂平均, n∈(0,∞)) ───
F_Ricci = ( f₁^n + f₂^n + ... + fₖ^n )^(1/n)    当所有 fᵢ > 0
F_Ricci = max(f₁,f₂,...,fₖ)                其他情况
// n→∞ → CSG 并 (max)；n→0 → 几何平均；n=1 → 算术平均

// ─── ③ 平滑最小 (Smooth Union / smin, 参数 k>0) ───
smin(f₁,f₂, k) = min(f₁,f₂) − max(k − |f₁−f₂|, 0)² / (4k)
// k 控制混合半径：|f₁−f₂|≥k 时退化为精确 min

// ─── ④ 平滑最大 (Smooth Intersection / smax) ───
smax(f₁,f₂, k) = max(f₁,f₂) + max(k − |f₁−f₂|, 0)² / (4k)

// ─── ⑤ R 函数 (α∈[−1,1], 保持精确布尔语义) ───
// R-并 (R-disjunction)：
f₁ ∨α f₂ = f₁ + f₂ + √( f₁² + f₂² − 2α·f₁·f₂ )
// R-交 (R-conjunction)：
f₁ ∧α f₂ = f₁ + f₂ − √( f₁² + f₂² − 2α·f₁·f₂ )
// α=1 退化为 2·max 和 2·min (CSG)

// ─── ⑥ Wyvill 可控混合 ───
F_Wyvill = f₁ + f₂ − √( f₁² + f₂² − k·f₁·f₂ )     k ∈ [0,2]

21.4.1 求和混合的涨大效应

求和混合 F = Σ f˰ − T 虽简单，但有严重缺陷——涨大效应（Bulging）。当两 blob 靠近时，重叠区场叠加使等值面向外膨胀，V形内角（如分枝处）尤甚。原因是线性叠加不能区分"两个物体各自独立"与"两个物体融合"的语义差异。

需要可控混合算子，在保持平滑连接的同时抑制非预期膨胀。

21.4.2 Ricci 混合

Ricci 混合（Ricci, 1973）通过幂平均参数化控制融合程度：

参数 n 的行为分析：

Ricci 混合的局限性：(1) 仅在 f˰>0（内部）混合，在外部（f˰<0）退化为 max，导致混合区域内外不对称；(2) 对三个及以上图元需限定全部为正的区域；(3) 不能指定混合区域的精确范围。

21.4.3 平滑 min 与平滑 max

直接从 CSG min/max 引入平滑过渡将产生更自然的混合：

这是基于二次多项式的 C¹ 平滑逼近，由 Quilez（iq）在 Shadertoy 社区广泛推广。其核心优势是局部混合控制——k 精确指定了被混合区域的宽度，超出此宽度完全退化为离散 CSG。

21.4.4 R 函数（R-functions）

R 函数（R-functions, Rvachev, 1963）是布尔逻辑的连续函数对应。若实函数 F(f₁, f₂) 的符号仅依赖于 f₁、f₂ 的符号（而非具体值），则 F 称为 R 函数。

退化情况验证：取 α=1，平方根内变为 (f₁−f₂)²：

R 函数保证：sign(F(f₁,f₂)) 完全由 sign(f₁) 和 sign(f₂) 的布尔逻辑决定。这是其区别于一般平滑算子的本质特征——保持精确的布尔语义，同时提供 C¹ 连续的几何过渡。

21.4.5 Wyvill 混合公式族

Wyvill 等人（1999）在 BlobTree 框架中开发了专用混合公式族：

21.4.6 混合算子速查对比表

21.5 精确多边形化

21.5.1 Marching Cubes 算法流程

渲染管线需要三角形。Marching Cubes（Lorensen & Cline, 1987）是最经典的隐式曲面多边形化算法，分四步：

体素顶点编号约定（局部坐标系，单位立方体 (0,0,0)-(1,1,1)）：

12 条边的端点定义（边 e 连接 topoEdge[e][0] 和 topoEdge[e][1]）：

══════════════════════════════════════════════
  Marching Cubes — 配置编码与边表结构
══════════════════════════════════════════════

// 8 顶点场值符号编码为 8 位无符号整数 (cube_index ∈ [0,255])
// 位序：bit(v) = f(v)<0 ? 1 : 0，v=0..7 按 Morton 序

cube_index = Σ_{v=0}^{7} ( f(vertex_v) < 0 ? 1 : 0 ) · 2^v

// 边表 edgeTable[256] 每项为 12 位位掩码：
//   bit(e) = 1 表示等值面与边 e 相交
//   格式：edgeTable[idx] = Σ_{e=0}^{11} (与边e相交 ? 1 : 0)·2^e

// 三角表 triTable[256] 每项为不定长边三元组列表：
//   triTable[idx] = {e_a, e_b, e_c,  e_d, e_e, e_f,  ... , -1}
//   每三个连续边索引定义一个三角形，-1 为终止哨兵

// 对称性简化：通过旋转 (4-fold)×反射 (2-fold) 共 8 种对称
//   将 256 种配置归约为 15 种基本拓扑 (Case 0-14)
//   实际实现可用 256 项全表 (无运行时对称变换，查表直出结果)

// 顶点法线计算 — 中心差分 (Central Difference)：
// ∇f(p) ≈ ( f(p+h·e_x)−f(p−h·e_x) ,
//            f(p+h·e_y)−f(p−h·e_y) ,
//            f(p+h·e_z)−f(p−h·e_z) ) / (2h)
// h = 体素边长 (或更小的偏移) 用于有限差分求梯度

21.5.2 15 种基本三角剖分配置（完整表）

8 个顶点的符号模式有 2⁸ = 256 种，通过旋转和反射对称性归约为 15 种基本配置。每种配置包含正顶点数、三角形数、涉及的边及三角形连接方式。

对称关系：Case 1 和 Case 14 互为补（正负反转），Case 2 和 Case 13（部分），Case 4 和 Case 11（部分）等。通过旋转 90°/180° 可生成每种配置的 4~24 种定向变体，对应全部 256 个 cube_index。

21.5.3 边上的线性插值

══════════════════════════════════════════════
  Marching Cubes — 边零点插值 (Edge Interpolation)
══════════════════════════════════════════════

// 问题：边 e 两端点 v₀,v₁ 的场值为 f₀=f(v₀) < 0, f₁=f(v₁) > 0
// 目标：找到边上 f(p)=0 的精确位置 (零点)
// 假定：标量场 f 在边上线性变化 (当体素足够小时这是良好近似)

// 线性插值：f(t) = (1−t)·f₀ + t·f₁
// 令 f(t)=0 求解 t：

           f₀
t = ───────────
       f₀ − f₁

// 交点位置 (世界坐标)：
p = (1 − t) · v₀ + t · v₁

// 推导验证：
//   f₀<0 → 分子为负 → t>0 ✓ (交点在 v₀ 右侧)
//   f₁>0 → 分母 f₀−f₁<0 → t<1 ✓ (交点在 v₁ 左侧)

//   t=0 时 p=v₀ (f₀=0 的退化情况，交点在端点)
//   t=1 时 p=v₁ (f₁=0 的退化情况，交点在端点)

// 数值稳定性：当 |f₀−f₁| < ε 时，边平行于等值面 → 跳过该边

设边 e 的两端点 v₀、v₁ 的场值为 f₀=f(v₀)、f₁=f(v₁)，其中 f₀ < 0、f₁ > 0（或反之），则零点位置插值：

推导：线性插值 f(t) = (1−t)f₀ + t·f₁，令 f(t)=0 解得 t=-f₀/(f₁−f₀) = f₀/(f₀−f₁)。

21.5.4 面歧义与内部体歧义

面歧义（Face Ambiguity）是 Marching Cubes 的核心难题。当体素某面的四个顶点中两个对角为正、两个对角为负时，该面上等值线有两种合法连接方式：

渐近判定器（Asymptotic Decider, Nielson & Hamann, 1991）通过检查面中心点的场值确定正确拓扑：

面中心点处的场值通过对角线双线性插值计算：

内部体歧义（Internal Ambiguity）出现在 Case 4/7/8/10/12/13 等配置中：即使面歧义已解决，体素内部仍可能有两种合法三角剖分（连通 vs 不连通的子分量）。需要额外采样体素中心点并比较场值符号来判定。

21.5.5 自适应 Marching Cubes

均匀分辨率有两难：平坦区浪费三角形，高曲率区分辨率不足。

══════════════════════════════════════════════
  自适应 Marching Cubes — 细化准则与判定
══════════════════════════════════════════════

// 细化判定准则 1：拓扑变化 (体素跨越零等值面)
//   若体素 8 顶点场值不全同号 → 等值面穿过 → 需要细化
判定: min(f₀..f₇) · max(f₀..f₇) < 0

// 细化判定准则 2：曲率自适应 (高曲率处提高分辨率)
//   计算体素中心点的平均曲率估计：
κ ≈ ||∇²f|| / ||∇f||
//   若 κ > κ_max / voxel_size → 继续细化 (高曲率需更细体素)

// 细化判定准则 3：特征大小 (防子体素细节丢失)
//   在体素内部随机采样 N 点，检查是否 f 穿过零：
判定: f(pₛ) 的符号在采样点间不一致 → 可能丢失细节 → 细化

// ─── T-junction 修补：限制八叉树 (Restricted Octree) ───
// 约束：相邻叶节点层级差 ≤ 1 (2:1 平衡条件)
// 此时 T-junction 退化为有限情况，可用三角剖分补丁处理

// ─── 自适应 vs 均匀网格 — 三角形数对比 ───
//   均匀 128³ 网格 → ~2.1M 体素 → ~500K 三角形 (含大量平坦区浪费)
//   自适应至 8 层 (最小 1/256) → ~50K 体素 → ~120K 三角形 (省 75%)

自适应 Marching Cubes（Adaptive Marching Cubes, AMC）基于八叉树：在零等值面曲率大处细分体素、平坦处保持粗粒度。关键挑战：不同分辨率邻居间的 T-junction 裂缝。

T-junction 的产生与修补

当一个体素边被细分而相邻体素的同一边未细分时，两三角形间产生 T 形接缝——若不加处理，该接缝处顶点只属于一侧三角形而不属于另一侧，导致可视化裂缝。

21.5.6 其他等值面提取方法

Marching Cubes 并非唯一选择。几种演进算法针对其弱点做了改进：

在 BlobTree 框架下，Marching Cubes 仍是默认选择——实现成熟、查表完备、与三角形渲染管线直接对接。Dual Contouring 在需要保留由 CSG 操作产生的尖锐边角时表现更好。

21.6 变形

21.6.1 空间变形的概念

隐式建模一大优势：对空间本身施加变形，而非逐个顶点计算。设变形映射 D : ℝ³ → ℝ³：

重要：这里 D 是"逆映射"——对 p 先做 D(p) 得到原空间坐标，再用原隐式函数求值。等价于原物体被逆映射 D−¹ = {q | D(q)=p} 变形到新位置。

变形的核心优势——与几何复杂度解耦：无论物体有多少细节，变形操作仅在对点 p 求值时额外计算一次 D(p) 及其导数，代价 O(1)。这是隐式建模对有机形状建模最有力的武器。

21.6.2 常见全局变形公式

Twist（扭转）

绕 z 轴按高度比例扭转——离原点越远扭转角度越大：

雅可比矩阵（用于法线修正）：

Bend（弯曲）

沿 x 轴将空间弯曲成弧形（绕 y 轴弯曲）：

检查一致性：在 x=0 处，sin(0)=0，cos(0)=1，故 D(0,y,z) = (0, y, z+0) = (0,y,z)——不变。小 δ 偏移 x 对应弧长 δ——线性近似成立。

雅可比矩阵：

det(J_D) = cos(x/R) > 0 当 |x| < πR/2——变形在有效弯曲角度内保持定向。

Taper（锥化）

沿 z 轴线性缩放 xy 分量，使物体一端粗一端细：

Shear（剪切）

沿某个方向按另一方向的坐标偏移：

组合变形

多个变形可通过函数复合实现复杂效果：

21.6.3 全局变形 vs 局部变形

局部变形需额外的衰减函数（Falloff Function）调节变形强度。设影响区域 Ω：

常用衰减函数：

21.6.4 变形后的法线修正

══════════════════════════════════════════════
  空间变形 — 雅可比矩阵与法线修正
══════════════════════════════════════════════

// 变形映射 D : ℝ³ → ℝ³ 定义逆变形求值：
F'(p) = F( D(p) )

// ─── Twist (绕 z 轴扭转，θ=α·z) ───
// D(p) = ( x·cosθ + y·sinθ , −x·sinθ + y·cosθ , z )
// 雅可比 J_D = ∂D/∂p：
       ┌                                      ┐
       │  cosθ        sinθ       −α(x·sinθ−y·cosθ) │
J_D = │ −sinθ       cosθ       −α(x·cosθ+y·sinθ) │
       │   0           0               1           │
       └                                      ┘
// det(J_D) = cos²θ + sin²θ = 1 → 体积保持 (等容变形)

// ─── Bend (沿 x 轴弯曲，弯曲半径 R) ───
// D(p) = ( R·sin(x/R) , y , R−R·cos(x/R)+z )
       ┌                       ┐
       │  cos(x/R)   0   0     │
J_D = │    0        1   0     │
       │  sin(x/R)   0   1     │
       └                       ┘
// det(J_D) = cos(x/R) > 0 当 |x|<πR/2 (保持定向)

// ─── Taper (沿 z 轴锥化，s(z)=1−k·z) ───
// D(p) = ( x/s(z) , y/s(z) , z )
       ┌                                    ┐
       │  1/s(z)    0      k·x/s(z)²       │
J_D = │    0     1/s(z)   k·y/s(z)²       │
       │    0       0         1            │
       └                                    ┘

// ─── 变形后法线修正 (Normal Correction) ───
// 核心公式：法线按逆雅可比的转置变换

n'(p) = ( J_D⁻¹(p) )ᵀ · ∇F( D(p) )

// 推导：由链式法则
//   ∇F'(p) = J_D(p)ᵀ · ∇F(D(p))  →  n'(p) = (J_D⁻¹)ᵀ · n(D(p))
// 原因：法线是余向量 (covector)，按逆转置变换以保持垂直于切平面

// ─── 组合变形 (复合顺序重要) ───
F'(p) = F( D₁( D₂( ... Dₖ(p) ... ) ) )
// 对应雅可比：J = J_D₁ · J_D₂ · ... · J_Dₖ (矩阵乘法，不可交换)
// 先 Twist 再 Bend ≠ 先 Bend 再 Twist

变形 D 不仅移动空间位置，还扭曲局部几何框架。原物体的法向量 n 需按 D 的逆雅可比转置修正：

推导原理：由链式法则，∇F'(p) = J_D−¹(p)ᵁ · ∇F(D(p))。在物体曲面保持局部体积时，法线变换需用逆转置保持垂直于切平面。

21.7 参数化

══════════════════════════════════════════════
  隐式曲面参数化 — 公式与方法汇编
══════════════════════════════════════════════

// 问题：给定隐式曲面 S = {p | f(p)=0}，求映射 φ: S → Ω ⊂ ℝ²
// 困难：隐式表示无内建参数坐标 (与 p(u,v) 参数曲面不同)

// ─── 方案 1：投影参数化 ───
// 平面投影 (基向量 e₁, e₂)：
   (u, v) = ( p·e₁ , p·e₂ )
// 柱面投影 (沿 y 轴)：
   (u, v) = ( atan2(p_z, p_x) , p_y )
// 球面投影：
   (u, v) = ( atan2(p_z, p_x)/2π , acos(p_y/||p||)/π )

// ─── 方案 2：LSCM (最小二乘保角映射, Lévy et al. 2002) ───
// 在已生成的三角网格上求解保角参数化
// 目标函数：每个三角形 T 上的梯度应满足 Cauchy-Riemann 关系

   min_{u,v}  Σ_{T} || ∇u(T) − i · ∇v(T) ||² · A(T)

// 展开为离散形式 (三角形 T 顶点 1,2,3)：
//   ∂u/∂x ≈ (u₂−u₁)(y₃−y₁)−(u₃−u₁)(y₂−y₁) / 2A(T)
//   ∂u/∂y ≈ (u₃−u₁)(x₂−x₁)−(u₂−u₁)(x₃−x₁) / 2A(T)
// 形成稀疏线性系统 Au = 0 (需固定至少 2 个顶点消歧义)

// ─── 方案 3：骨架属性传递 (BlobTree 原生方案) ───
// 每叶子图元储存其局部 (u,v) 映射，混合节点按场值加权：

   attr(p) = Σᵢ wᵢ(p) · attrᵢ(p) / Σᵢ wᵢ(p)
   wᵢ(p) = max( fᵢ(p) , 0 )

// ─── 方案 4：活动标架法 (Moving Frame) ───
// 基于梯度场构建局部切向量基，积分得到 (u,v)
   n̂ = ∇f/||∇f||                     法向量
   û = (n̂ × e_z)/||...||              第一切向量 (取与 z 轴叉积)
   v̂ = n̂ × û                       第二切向量 (右手法则)
// 受毛球定理制约：封闭曲面至少 2 个奇点 → 需割缝

21.7.1 隐式到参数化：根本挑战

纹理映射、UV 展开、烘焙等需要参数域 (u,v) 的二维坐标。从隐式曲面（定义为由 f(p)=0 的隐式点集）到参数化映射的过渡面临本质困难——隐式表示不提供内在的参数坐标。

形式上，需要找到映射 φ : S → Ω ⊂ ℝ²，其中 S = {p | f(p)=0}。与参数曲面不同（参数域天然由 u,v 定义域给出），隐式曲面需额外构造该映射——这就是隐式→参数化反问题（Inverse Problem）。

21.7.2 实用方案

方案 1：投影参数化

将曲面投影到简单形状上生成 (u,v)：

投影法快速但不保角也不保面积——在曲面法线与投影方向接近垂直处产生奇点（纹理拉伸/挤压）。

方案 2：离散参数化

先 Marching Cubes 获得三角形网格 M ≈ S，再应用标准网格参数化算法：

优缺点：网格质量依赖分辨率，但可获得高质量 UV；额外的网格生成开销。

方案 3：骨架属性传递

在 BlobTree 框架下，参数附着在骨架图元上。每个图元储存其局部 (u,v) 映射（对球体用经纬坐标，圆柱体用轴向+周向坐标等）。混合节点定义组合规则——按场值权重混合 (u,v)：

这种方法的优势是参数化与场求值并行——无需后处理步骤，且自然地处理了混合区域的参数过渡。

21.7.3 梯度方向场与活动标架

隐式场梯度 ∇F(p) 给出曲面法向量和局部方向信息。可由此构建曲面上的活动标架（Moving Frame）以定义局部 (u,v)：

但这一局部标架受毛球定理（Hairy Ball Theorem）制约——曲面上不可能存在处处非零的连续切向量场。对亏格为 0 的封闭曲面（拓扑球面），任何切向量场至少有两个零点。这意味着基于梯度的全局参数化必然有奇点，需通过割缝（Seam）或图集（Atlas）来解决。

21.8 BlobTree

21.8.1 层次化隐式场景结构

BlobTree（Wyvill et al., 1999）将隐式建模的所有操作统一到树形数据结构中。叶子节点是骨架图元（生成基础标量场），内部节点是操作（混合、CSG、变形），根节点输出最终标量场 F(p)。

21.8.2 完整节点类型

21.8.3 递归求值遍历

BlobTree 的求值是自底向上的递归遍历：

加速策略：空间划分结构（21.3 节）在求值前剪枝——仅对支撑域覆盖 p 的叶子图元求值。对混合节点，若某子节点返回 0 或负值可提前终止（对 min 操作）。

══════════════════════════════════════════════
  BlobTree — 节点求值数学定义
══════════════════════════════════════════════

// 每种节点类型的精确求值公式：

// Primitive (叶子节点)：直接计算 Wyvill 场
f_prim(p) = (1 − ||p−c||²/R²)³     当 ||p−c|| ≤ R，否则 0

// Union (CSG 并)：场值取最大 (内部区域取并)
f_∪(p) = max{ fᵢ(p) }                     子节点 2..k 个

// Intersection (CSG 交)：场值取最小 (内部区域取交)
f_∩(p) = min{ fᵢ(p) }                     子节点 2..k 个

// Difference (CSG 差 A−B)：A 中挖去 B 的内部区域
f_−(p) = min( f_A(p), −f_B(p) )

// Deform (空间变形)：对查询点 p 先施加逆变形再求子节点场
f_deform(p) = child( D(p) )              D 是逆变形映射

// Blend (混合节点, 根据类型选择公式)：
//   Sum:    F = Σ fᵢ − T
//   Ricci:  F = (Σ fᵢ^n)^(1/n)  (当 fᵢ>0，否则 max)
//   Smooth: F = smin/smax(f₁,f₂,k)
//   Wyvill: F = f₁+f₂−√(f₁²+f₂²−k·f₁·f₂)

// 树求值复杂度 (含空间加速)：
//   无加速：O(L) 每点求值 (L=叶节点数)
//   八叉树：O(K·log N) 每点求值 (K=每单元平均图元数)

21.8.4 射线追踪渲染隐式曲面：光线步进算法

BlobTree 可通过光线步进（Ray Marching）直接渲染，无需先生成三角形网格。这是隐式建模最"原生"的渲染方式。

SDF 安全步长原理：当使用 SDF 时，|f(p)| 即为 p 到曲面的最短距离。因此可以安全地沿光线方向前进 |f(p)| 而不会越过曲面（球追踪 Sphere Tracing, Hart 1996）。非 SDF 场（如 Wyvill 求和场）不具备单位梯度性质，只能使用固定步长或自适应步长，效率较低。

══════════════════════════════════════════════
  球追踪 (Sphere Tracing) — 数学原理与步长分析
══════════════════════════════════════════════

// 球追踪的核心：对 SDF，|f(p)| = dist(p, ∂Ω)
// 沿光线方向前进 |f(p)| 保证不越过曲面

// 证明：设 q 为 ∂Ω 上距 p 最近的点
//   则球 B(p, |f(p)|) 内部无曲面点
//   光线在该球内任意前进 ≤ |f(p)| 都不会穿过 ∂Ω
//   在 p + step·dir 处，f 的符号尚未改变 → 仍在同侧

// 步长策略对比：
//   SDF 球追踪：  t += |f(p)|              → 每步自适应
//   固定步长：    t += Δt                  → Δt=0.001~0.01
//   Lipschitz：   t += |f(p)|/L            → L=sup ||∇f||
//   过松弛球追踪：t += ω·|f(p)|  (ω∈[1,2]) → 加速但有越过风险

// ─── 光线步进总步数分析 (从 d₀=10 到球面 d=1) ───
//   SDF 球追踪：每步至少前进 |f(p)|，步数 ≈ log(d₀/ε) ≈ 10-15
//   固定步长 0.01：步数 = 9/0.01 = 900
//   Lipschitz (Wyvill L≈3/R)：步数 ≈ 3·log(d₀/ε) ≈ 30-45
//   → SDF 球追踪比固定步长快 60-90×！

// ─── 二分法求精 (Bisection Refinement) ───
// 光线步进发现 f(p_prev)>0 and f(p_curr)<0 后：
   binary_search(p_prev, p_curr, ε):
     for i = 0..10:
       t_mid = (t_prev + t_curr) / 2
       if |f(p_mid)| < ε: return t_mid
       if f(p_mid) < 0: t_curr = t_mid
       else: t_prev = t_mid
     return (t_prev + t_curr) / 2
// 10 次二分将误差从体素边长缩至 体素边长/1024 ≈ 0.001

Lipschitz 界与自适应步长

对于非 SDF 场 f，可估计其Lipschitz 常数 L（梯度的模上界）：

保守步长可取 |f(p)| / L，比固定步长更安全但可能过于保守（因 L 是全局上界，局部 |∇f| 可能远小于 L）。

21.8.5 属性传递：颜色、材质与更多

BlobTree 可传递除标量场外的任意属性——颜色、法线、材质参数等。每个叶子图元带有其局部属性（球体天然经纬色、圆柱轴向色等），内部节点定义组合规则。

对 CSG 节点，被裁剪区域（非正场子节点）不贡献属性，这自然符合布尔操作的语义。

21.9 交互式隐式建模

══════════════════════════════════════════════
  交互式隐式建模 — 性能模型与加速公式
══════════════════════════════════════════════

// ─── 实时预览性能模型 ───
// 交互延迟 = 场求值延迟 + 曲面提取延迟 + 渲染延迟
// 场求值延迟 (N=场景总图元数, K=活跃单元平均图元数)：
   T_eval = O(K) · cost(Wyvill)        含八叉树加速
// 曲面提取延迟 (M=体素总数, G=跨零体素数)：
   T_extract = M · O(1) + G · O(1)    每体素 8 次求值 + 跨零体素查表

// ─── Wyvill 查表加速 (LUT) ───
// 预计算 (1−t²)³ 在 t∈[0,1] 上 256 等分采样：
   table[k] = (1 − (k/255)²)³      k = 0,1,...,255
// 运行时：t = r/R → idx = floor(t·255) → 查表 + 线性插值
//   多项式求值 (~12 FLOP) → 查表+LERP (~3 FLOP) → ~4× 加速

// ─── GPU 并行 Marching Cubes ───
// Compute Shader 中每线程处理一个体素：
//   f[0..7] ← Wyvill(vertex_0)..Wyvill(vertex_7)   // 8 次求值
//   cube_idx ← pack_signs(f[0..7])                   // 压缩为 8 位
//   edges ← edgeTable[cube_idx]                      // 查表
//   for each tri in triTable[cube_idx]: emit_tri()   // 写三角形
//   三角形写入 AppendStructuredBuffer (原子计数分配)

// ─── Lipschitz 自适应步长 ───
// 对非 SDF 隐式场，估计最大梯度上界 L：
   L = sup_{p∈Ω} ||∇f(p)||
// 保守步长：step ≤ |f(p)| / L
// 对 Wyvill 函数：L ≈ 3/R (在 r=R/√3 处梯度最大)

21.9.1 实时预览策略

隐式建模的交互瓶颈在于：每次参数修改需要重新求值全场（或受影响子区域）并重建曲面（Marching Cubes 或光线步进）。策略：

增量刷新区

仅重算受影响的局部区域（脏矩形 Dirty Region）：

渐进式细化

先低分辨率预览，后台逐步提高分辨率。用户拖动图元时用 32³ 体素网格粗预览，释放鼠标后逐步细化至 128³ 甚至 256³。多数现代隐式建模工具都采用此策略。

GPU 加速

GPU Marching Cubes：在计算着色器中批量求值体素顶点并生成三角形。
GPU 光线步进：片段着色器中直接求值 BlobTree，像素级精度渲染——无需任何三角形化。
混合方案：GPU 计算着色器做 Marching Cubes + 光栅化渲染最终三角形。

21.9.2 雕刻系统隐式工具

交互式隐式建模中的雕刻工具（Sculpting Tools）使艺术家像处理黏土一样工作：

工具 1：添加 blob（Add）

点击/拖拽在空间中添加新的骨架图元，自动与现有场混合：

工具 2：减去 blob（Subtract）

添加负场 blob 实现"挖洞"效果：

工具 3：平滑（Smooth）

对局部区域降低场梯度，平滑曲面：

工具 4：局部变形笔刷（Deform Brush）

在影响区域内拉动/推压曲面点：

工具 5：CSG 剪辑器

用隐式几何（如球、长方体、圆柱）进行布尔并/交/差裁剪——结果即时可见。因 CSG 操作在 O(1) 内完成（单次求值中的 min/max），这是隐式建模在概念建模阶段超越多边形布尔的最大优势。

21.9.3 用户界面范式

21.9.4 性能优化要点

课后练习题（7 题）

1. Wyvill 导数与连续性：证明 f(r)=(1−r²)³ 在 r=1 处 C² 连续（f、f'、f'' 均连续且为零）。求 r=0.5 处的场值和梯度大小。若半径 R=2，计算 r=1.5 处的场值。
2. Blinn vs Wyvill：对于 Blinn blob f(p)=exp(−2·r²)−0.1，有效半径 R 为多少？在 r=R 处梯度模 ||∇f|| 为多少？与 Wyvill 函数比较，两种模型在求值效率和混合行为上有何差异？
3. Matching Cubes 歧义排查：取体素顶点 (0,0,0) 处场值为 +1，(0,1,1) 处为 +1，(1,0,1) 为 −1，(1,1,0) 为 −1（四个正、四个负，呈四面体排列）。判断此配置对应基本配置表哪个 case。相邻体素不同剖分的后果是什么？
4. Ricci 混合极限分析：两相等球场 f₁=f₂=1−r²，中心间距 d=√2。令 Ricci 混合 F=(f₁ᵗ+f₂ᵗ)¹ˆᵗ，画 n=1,2,5,50 时两中心连心线上的 F 值剖面图，分析膨胀程度随 n 的变化。当 n→∞ 时零等值面何去何从？
5. R 函数验证：验证 R 交 ∧ᵦ 在 α=1 时退化为 2·min(f₁, f₂)，证明零等值面对应于 f₁≤0 且 f₂≤0 区域的边界——正确实现了布尔交的语义。
6. BlobTree 复杂度分析：深度 h、叶节点 L 的 BlobTree。每变形节点 O(1)（矩阵乘法），混合节点 m 子节点。无空间加速时推导单点求值最坏时间复杂度。加入八叉树后复杂度如何变化？
7. 光线步进步长策略对比：对 SDF 球体 f(p)=||p||−1，从 (0,0,−10) 沿 (0,0,1) 方向：计算球追踪（安全步长）和固定步长（Δt=0.1）所需的步数。扩展到一般隐式场（非 SDF）讨论自适应步长的 Lipschitz 界。

QA 零基础问答区（8 问）

— FCG 第5版 · 第21章：隐式建模 · 全文翻译（充实的级别） —

Guizang 插图：Marching Cubes 配置表 · Wyvill BlobTree · 变形公式 © 相关版权归原作者所有 · 讲义仅供学习使用