Luna Flow
简体中文

core 设计

设计目标

core 是 geometry3d 中所有其他包都认同的部分：网格是什么，点在变换下如何移动，面的哪一侧是外侧，以及面在光照下有多亮。它必须小到一次就能读完，对约定准确无误，并且不包含任何属于相机、屏幕或输出设备的东西。它也展示了 Luna-Flow/linear-algebra 的稠密类型足以支撑一条 3D 管线。

数学背景

齐次坐标

R3\mathbb{R}^3 上的仿射映射 x↦Ax+tx \mapsto A x + t 不是线性的，因此没有 3×3 矩阵。把 R3\mathbb{R}^3 嵌入 R4\mathbb{R}^4 就能解决这个问题。点 pp 写作 (p,1)(p, 1)，方向 dd 写作 (d,0)(d, 0)，于是仿射映射变成分块矩阵

M=(At0T1),M(p1)=(Ap+t1),M(d0)=(Ad0).M = \begin{pmatrix} A & t \\ 0^\mathsf{T} & 1 \end{pmatrix},\qquad M \begin{pmatrix} p \\ 1 \end{pmatrix} = \begin{pmatrix} A p + t \\ 1 \end{pmatrix},\qquad M \begin{pmatrix} d \\ 0 \end{pmatrix} = \begin{pmatrix} A d \\ 0 \end{pmatrix}.

最后一个坐标记录了向量是哪种对象。方向是两点之差，(p,1)−(q,1)=(p−q,0)(p, 1) - (q, 1) = (p - q, 0)，而上面的计算表明它只受 AA 作用：平移两个点不会改变它们的差。这正是 Transform3::apply_point（追加 w=1w = 1）与 Transform3::apply_direction（追加 w=0w = 0）之间的区别。

一般的 4×4 矩阵还可能有最后一行 (hT,k)≠(0,0,0,1)(h^\mathsf{T}, k) \ne (0, 0, 0, 1)。此时它把 (p,1)(p, 1) 送到 (x′,y′,z′,w′)(x', y', z', w')，其中 w′=h⋅p+kw' = h \cdot p + k，而它所表示的点通过齐次除法 (x′/w′,y′/w′,z′/w′)(x'/w', y'/w', z'/w') 得到。由于对每个 λ≠0\lambda \ne 0，(x′,y′,z′,w′)(x', y', z', w') 与 λ(x′,y′,z′,w′)\lambda (x', y', z', w') 给出同一个点，这些矩阵作用在射影空间上；其中包括 view 所用的透视投影。apply_point 总是做除法，因此两种情况都能处理，并在 ∣w′∣≤ε|w'| \le \varepsilon = DEPTH_EPSILON 时，即对被送到（或接近）无穷远平面的点，跳过除法。

组合

对于两个变换 M1M_1（先应用）和 M2M_2（后应用），由矩阵乘法的结合律，复合变换作用于列向量为 M2(M1v)=(M2M1)vM_2 (M_1 v) = (M_2 M_1) v。Transform3::compose 存储 M2M1M_2 M_1：

t1.compose(t2).matrix=M2M1.\texttt{t1.compose(t2)}.\mathrm{matrix} = M_2 M_1 .

该方法按应用顺序阅读，而矩阵乘积从右往左阅读。由代数可得两个结论。组合满足结合律，因此链可以任意分组。它不满足交换律：设 TT 为平移 tt、RR 为旋转，

TR(p1)=(Rp+t1)≠RT(p1)=(Rp+Rt1)T R \begin{pmatrix} p \\ 1 \end{pmatrix} = \begin{pmatrix} R p + t \\ 1 \end{pmatrix} \quad\ne\quad R T \begin{pmatrix} p \\ 1 \end{pmatrix} = \begin{pmatrix} R p + R t \\ 1 \end{pmatrix}

除非 Rt=tR t = t。因此通常的模型变换先缩放，再旋转，最后平移：s.compose(r).compose(t)，矩阵为 TRST R S。

基本旋转

线性映射由基向量的像决定，这些像构成其矩阵的各列。把 xyxy 平面绕 zz 轴旋转 θ\theta，会把

ex↦(cos⁡θ,sin⁡θ,0),ey↦(−sin⁡θ,cos⁡θ,0),ez↦ez,e_x \mapsto (\cos\theta, \sin\theta, 0),\qquad e_y \mapsto (-\sin\theta, \cos\theta, 0),\qquad e_z \mapsto e_z ,

由此得到 rotation_z 的矩阵 Rz(θ)R_z(\theta)。绕 xx 和 yy 的旋转可通过循环轮换坐标轴 x→y→z→xx \to y \to z \to x 得到：绕 xx 时由 (y,z)(y, z) 扮演 (x,y)(x, y) 的角色，绕 yy 时由 (z,x)(z, x) 扮演。第二次替换解释了为何 RyR_y 的负号位于对角线下方：

Ry(θ):ez↦(sin⁡θ,0,cos⁡θ),ex↦(cos⁡θ,0,−sin⁡θ).R_y(\theta):\quad e_z \mapsto (\sin\theta, 0, \cos\theta),\qquad e_x \mapsto (\cos\theta, 0, -\sin\theta).

每个 RR 都是行列式为 +1+1 的正交矩阵：其各列标准正交，且 det⁡Rz(θ)=cos⁡2θ+sin⁡2θ=1\det R_z(\theta) = \cos^2\theta + \sin^2\theta = 1。因此 R−1=RT=R(−θ)R^{-1} = R^\mathsf{T} = R(-\theta)，旋转保持长度、角度和标架的定向（把右手三元组映射为右手三元组）。

欧拉角

rotation_matrix(α, β, γ) 把基本旋转组合为 R=Rz(γ)Ry(β)Rx(α)R = R_z(\gamma) R_y(\beta) R_x(\alpha)。记 c∙=cos⁡c_\bullet = \cos、s∙=sin⁡s_\bullet = \sin（作用于对应角度）并展开相乘，

Ry(β)Rx(α)=(cβsβsαsβcα0cα−sα−sβcβsαcβcα),R=(cγcβcγsβsα−sγcαcγsβcα+sγsαsγcβsγsβsα+cγcαsγsβcα−cγsα−sβcβsαcβcα).\begin{aligned} R_y(\beta) R_x(\alpha) &= \begin{pmatrix} c_\beta & s_\beta s_\alpha & s_\beta c_\alpha \\ 0 & c_\alpha & -s_\alpha \\ -s_\beta & c_\beta s_\alpha & c_\beta c_\alpha \end{pmatrix},\\[4pt] R &= \begin{pmatrix} c_\gamma c_\beta & c_\gamma s_\beta s_\alpha - s_\gamma c_\alpha & c_\gamma s_\beta c_\alpha + s_\gamma s_\alpha \\ s_\gamma c_\beta & s_\gamma s_\beta s_\alpha + c_\gamma c_\alpha & s_\gamma s_\beta c_\alpha - c_\gamma s_\alpha \\ -s_\beta & c_\beta s_\alpha & c_\beta c_\alpha \end{pmatrix}. \end{aligned}

作用于列向量时，先绕 xx 旋转，然后 yy，再 zz，每次都绕固定的世界坐标轴（外旋 xx-yy-zz，等价于内旋 zz-y′y'-x′′x''）。旋转的乘积仍是旋转，因此由欧拉旋转定理，RR 是绕某个轴旋转某个角度 θ\theta 的单个旋转，且 tr⁡R=1+2cos⁡θ\operatorname{tr} R = 1 + 2\cos\theta。11 迹在基变换下不变，而在以旋转轴为第三个基向量的基下，矩阵为 Rz(θ)R_z(\theta)，其迹为 1+2cos⁡θ1 + 2\cos\theta。geometry3d 不提取轴和角；这个恒等式可用于在测试中检查组合后的旋转。

欧拉角存在奇点。当 β=π/2\beta = \pi/2 时，矩阵变为

R=(0sin⁡(α−γ)cos⁡(α−γ)0cos⁡(α−γ)−sin⁡(α−γ)−100),R = \begin{pmatrix} 0 & \sin(\alpha - \gamma) & \cos(\alpha - \gamma) \\ 0 & \cos(\alpha - \gamma) & -\sin(\alpha - \gamma) \\ -1 & 0 & 0 \end{pmatrix},

它只依赖于 α−γ\alpha - \gamma：把 α\alpha 和 γ\gamma 改变相同的量会得到同一个旋转。失去了一个自由度（万向节锁），而在 β=±π/2\beta = \pm\pi/2 附近，姿态的微小变化需要角度的大幅变化。演示程序以相互独立的恒定速率驱动三个角度，此时这一点无关紧要。

面、法向量与定向

当面 (a,b,c,d)(a, b, c, d) 的四个顶点共面时，该面位于一个平面内。其法向量由前三个顶点计算：

n=(vb−va)×(vc−va),n^=n/∥n∥.n = (v_b - v_a) \times (v_c - v_a),\qquad \hat n = n / \lVert n \rVert .

叉积是反对称的，因此交换两个顶点会使 nn 反向；顶点顺序就是定向。生成器列出每个面时都使 n^\hat n 指向实体外部。对于位于 z=−sz = -s 的立方体面 (0,3,2,1)(0, 3, 2, 1)：

v3−v0=(0,2s,0),v2−v0=(2s,2s,0),(0,2s,0)×(2s,2s,0)=(2s⋅0−0⋅2s, 0⋅2s−0⋅0, 0⋅2s−2s⋅2s)=(0,0,−4s2),\begin{aligned} v_3 - v_0 &= (0, 2s, 0), \qquad v_2 - v_0 = (2s, 2s, 0),\\ (0, 2s, 0) \times (2s, 2s, 0) &= (2s \cdot 0 - 0 \cdot 2s,\ 0 \cdot 2s - 0 \cdot 0,\ 0 \cdot 2s - 2s \cdot 2s) = (0, 0, -4s^2), \end{aligned}

它指向 −z-z，即远离立方体的方向。对于圆环 P(u,v)=((R+rcos⁡v)cos⁡u, rsin⁡v, (R+rcos⁡v)sin⁡u)P(u, v) = ((R + r\cos v)\cos u,\ r \sin v,\ (R + r\cos v)\sin u)，位于 (u,v)(u, v) 的面使用顶点 P(u,v)P(u, v)、P(u,v+Δv)P(u, v + \Delta v)、P(u+Δu,v+Δv)P(u + \Delta u, v + \Delta v)。一阶近似下两条边为 PvΔvP_v \Delta v 和 PvΔv+PuΔuP_v \Delta v + P_u \Delta u，因此

n≈(Pv×Pu) Δu Δv,n \approx (P_v \times P_u)\, \Delta u\, \Delta v ,

而在 u=v=0u = v = 0 处，Pv=(0,r,0)P_v = (0, r, 0) 与 Pu=(0,0,R+r)P_u = (0, 0, R + r) 给出 Pv×Pu=(r(R+r),0,0)P_v \times P_u = (r (R + r), 0, 0)：法向量背离旋转轴，指向管外。仓库中有一个测试检查这一点，以免背面剔除退化为显示内壁。

生成器输出的每个四边形都是平面的，因此其法向量是良定义的。对于球和圆环，参数 u,u′u, u' 之间的面在关于过 yy 轴、角度为 (u+u′)/2(u + u')/2 的平面的反射下对称，该反射交换 P(u,v)↔P(u′,v)P(u, v) \leftrightarrow P(u', v) 以及 P(u,v′)↔P(u′,v′)P(u, v') \leftrightarrow P(u', v')。线段 P(u,v)P(u′,v)P(u, v)P(u', v) 和 P(u,v′)P(u′,v′)P(u, v')P(u', v') 都垂直于该镜面，因而平行，而两条平行线段张成一个平面。立方体、圆柱侧面和退化四边形按构造就是平面的。

背面可见性

点 ee 能看到平面面片的正面，当且仅当它位于法向量所指向的开半空间中：

n^⋅(e−q)>0for a point q of the plane.\hat n \cdot (e - q) > 0 \quad\text{for a point } q \text{ of the plane.}

符号与 qq 的选择无关：对平面上的两点 q,q′q, q'，有 n^⋅(q−q′)=0\hat n \cdot (q - q') = 0，所以 n^⋅(e−q)=n^⋅(e−q′)\hat n \cdot (e - q) = \hat n \cdot (e - q')。face_is_visible 取 qq = face_center。对于封闭实体，背向眼睛的面会被同一实体的正面遮挡，因此剔除它永远不会去掉可见表面；这大约让光栅化器的工作量减半。它不解决不同正面之间的遮挡；那是深度缓冲的任务。

朗伯着色

面积为 AA、单位法向量为 n^\hat n 的小表面片，被沿单位方向 ℓ\ell 反向传播的平行光照射时，截获的光通量等于穿过垂直于光线、面积为 Acos⁡θA \cos\theta 的区域的通量，其中 cos⁡θ=n^⋅ℓ\cos\theta = \hat n \cdot \ell。理想漫反射（朗伯）表面向各个方向均匀反射光，因此其表观亮度与 cos⁡θ\cos\theta 成正比，当光源在表面背后时为零：

I=max⁡(0, n^⋅ℓ)∈[0,1].I = \max(0,\ \hat n \cdot \ell) \in [0, 1] .

face_intensity 返回的正是这个值，这就是 light 必须是指向光源的单位向量的原因。每个面一个值，就得到平面（分面）着色。

设计决策

以四边形为标准拓扑

问题：网格需要一种所有生成器、剔除测试和光栅化器都理解的面类型。可选方案有仅三角形、仅四边形或一般多边形。选择四边形是因为生成器的参数曲面（球、圆柱、圆环）是 (u,v)(u, v) 上的网格，其单元是四边形；一个四边形有一个法向量和一个着色值，使平面着色的面数减半。三角形以 d=ad = a 的退化四边形嵌入，因此棱锥、圆锥侧面或球的极冠都适用同一类型。triangulate_quad 只在光栅化器需要时才把四边形转成两个三角形，退化四边形的零面积第二个三角形会被面积测试丢弃。

点与方向作为不同的操作

Transform3 存储一个 4×4 矩阵，并提供 apply_point 和 apply_direction，而不是让调用者自己构造 ww。这把齐次约定封装在包内部，使“方向忽略平移”这一规则在调用处不可能出错。

重新计算法向量，而不是变换它们

在非均匀缩放下，法向量的变换方式与方向不同：法向量的正确矩阵是 (A−1)T(A^{-1})^\mathsf{T}，因为曲面的切向量 tt 满足 n⋅t=0n \cdot t = 0，而在映射 t′=Att' = A t 之后，向量 n′=A−Tnn' = A^{-\mathsf{T}} n 仍满足 n′⋅t′=nTA−1At=0n' \cdot t' = n^\mathsf{T} A^{-1} A t = 0。core 不存储法向量再用该矩阵变换，而是每次都由变换后的顶点重新计算 n^\hat n（在 apply_mesh 之后调用 face_normal）。代价是每个面一次叉积；好处是对任何 det⁡A>0\det A > 0 的仿射映射（包括非均匀缩放）法向量总是正确的。det⁡A<0\det A < 0 的映射（反射）会反转环绕方向，使每个面内外颠倒。

选择欧拉角，而不是轴角或四元数

演示程序只需要让物体以相互独立的速率绕三个轴旋转，这正是欧拉角所表达的，而且每个矩阵只有三行。轴角构造（罗德里格斯公式）和四元数组合没有实现；它们属于 Luna-Flow/quaternion。在别处生成的任何旋转矩阵都可以通过 Transform3::from_matrix 传入。

共享的容差

DEPTH_EPSILON = 10−910^{-9} 用于管线中每个“这是零吗？”的判断：归一化、齐次除法、三角形面积，以及严格的深度测试 depth + ε < stored。单一常量使退化边界处的行为在各包之间保持一致。它是绝对容差，因此假定场景坐标的量级在 11 到 10310^3 之间；演示程序使用的单位约为 1。

正确性与不变量

  • 对仿射变换调用 apply_point 在浮点舍入意义下是精确的：w′=1w' = 1，因此不会引入除法误差。
  • 对仿射变换，compose 满足 t1.compose(t2).apply_point(p)=t2.apply_point(t1.apply_point(p))\texttt{t1.compose(t2).apply\_point}(p) = \texttt{t2.apply\_point}(\texttt{t1.apply\_point}(p))；对射影变换则在齐次缩放意义下成立。有一个测试用先平移后缩放检查了顺序。
  • 旋转矩阵是行列式为 11 的正交矩阵；它们的乘积仍是旋转。
  • 生成的网格是封闭的、以原点为中心的，面是平面的且法向量朝外；face_is_visible 和 face_intensity 依赖这一点。
  • normalize_vec 和 face_normal 从不除以小于 DEPTH_EPSILON 的数；退化输入得到零向量，使该面不可见也不受光。
  • 每个操作都分配新的向量或网格；没有函数会修改其参数。apply_mesh 与其输入共享面数组。

除生成器和 apply_mesh（与顶点数和面数成线性关系）外，每个函数都是常数时间。

被否决的方案

  • 专用的小向量类型（带 x、y、z 字段的 Vec3 结构体）会更快，并且在维度上类型安全，但会重复 linear-algebra 的功能。本仓库的意义在于构建在 Luna-Flow 基础之上，因此向量是 @la.Vector[Double]，维度只是一种约定。
  • 泛型标量（对任意 Field 的 Mesh[T]）被否决：三角函数、平方根和容差都假定为 Double，而且没有后端能使用其他类型。
  • 在网格中存储法向量 因上述原因被否决：重新计算代价很小，并且在变换之后总是正确的。
  • 以三角形网格为主要类型 会使每个平面着色的网格面变成两个法向量相同的面。

边界

core 不会：

  • 了解相机、投影、视口、终端、颜色或 DOM；
  • 允许包外代码构造任意网格：Mesh 和 QuadFace 的字段在包外是只读的，因此生成器是网格的唯一来源；
  • 加载或保存网格，或提供场景图、材质、纹理、物理或空间索引；
  • 计算平滑（逐顶点）法向量、裁剪几何体或进行相交测试；
  • 提供逆变换、轴角或四元数旋转。

Footnotes

  1. 迹在基变换下不变，而在以旋转轴为第三个基向量的基下，矩阵为 Rz(θ)R_z(\theta)，其迹为 1+2cos⁡θ1 + 2\cos\theta。geometry3d 不提取轴和角；这个恒等式可用于在测试中检查组合后的旋转。 ↩