3.1 模型矩阵
模型矩阵是变换管线的第一段:它把一个物体从自己的局部坐标系搬到世界坐标系。本卷把整条链拆成「模型 → 视图 → 投影」三段 4x4 矩阵,逐段给出每一项的来历与定点化方式。如果你更关心 OpenGL 侧把矩阵通过 uniform 喂给哪一级可编程阶段、以及 GL 核心模式的 API 契约,那属于 OpenGL 核心模式渲染管线 的范畴;本节关心的是同一批 4x4 数字在 CPU 上如何被整数地构造、相乘、验证,以及为什么「顺序错一步,整帧像素就变样」。
这份关心不是修辞。本卷的核心主张是一份场景规格、三种宿主消费、两条一致性纪律:同一份 config/scene_spec.json 在 Python 与 TypeScript 下产出逐字节相同的像素缓冲。而像素缓冲逐字节相同的前提,是顶点变换链上每一次乘法的舍入方式都相同。模型矩阵正是这条链的第一个环节,它的复合顺序与舍入语义一旦两边不一致,后面所有像素都会分叉。
三种基本变换的齐次形式
三维里的平移不是线性变换,没法写成 3x3 矩阵。把它塞进 4x4 的办法是齐次坐标:点写成 (x, y, z, 1),向量写成 (x, y, z, 0)。第 4 列的三个分量于是承载平移,而第 4 行的 (0, 0, 0, 1) 保证齐次分量恒为 1。三种基本变换在 Q16.16 下分别是:
平移 T(t): 缩放 S(s): 绕 Y 轴旋转 R(c, s):
[1 0 0 tx] [sx 0 0 0] [ c 0 s 0]
[0 1 0 ty] [0 sy 0 0] [ 0 1 0 0]
[0 0 1 tz] [0 0 sz 0] [-s 0 c 0]
[0 0 0 1 ] [0 0 0 1] [ 0 0 0 1]
这里的 1 是 ONE = 65536,tx/ty/tz 与 sx/sy/sz 都是 Q16.16 整数。旋转矩阵用整数 c = cosθ、s = sinθ 常量构造;本节的算例取 θ = 90°,即 c = 0、s = ONE,于是每一项都是精确整数。任意角的 cos/sin 定点化留到第 4 章。
复合顺序为什么不可交换
矩阵乘法不满足交换律:T·R ≠ R·T。这一点在图形学里不是抽象代数趣闻,而是「先缩放还是先旋转」会给出不同模型的直接后果。原因在于矩阵作用在右侧:(T·R·S)·p 的求值顺序是「先 S、再 R、最后 T」。
T·R:先旋转、再平移——物体在原地转,然后被搬到新位置。R·T:先平移、再旋转——物体先被搬走,然后连同它的位移一起绕原点转。
同一个点 p = (1, 2, 3),两种顺序给出不同的落点。下面的脚本把它算出来,同时也验证了「复合矩阵」与「逐点分步应用」在本算例下逐项精确一致。
定点实现:先缩放、再旋转、再平移
脚本只复用 sim/cglib.py 的 mat_mul / mat_vec4,三个基本变换矩阵在本脚本内构造。它一次跑完四件事:复合矩阵作用、逐点分步作用、不可交换对照、矩阵乘中间量的量级统计。
#!/usr/bin/env python3
# file: lantern/sim/cg_raster_03_1_model_matrix.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_raster_03_1_model_matrix.py
"""3.1 模型矩阵:平移/旋转/缩放的 4x4 齐次形式、复合顺序、定点量级。
全整数(Q16.16,Python int);复用 sim/cglib.py 的 mat_mul / mat_vec4。
旋转矩阵在本脚本内用整数 cos/sin 常量构造(任意角三角函数留到第 4 章)。
"""
import math
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ONE = cglib.ONE
def mat_translate(t):
"""平移:平移量直接落在第 4 列(已是 Q16.16)。"""
return (ONE, 0, 0, t[0], 0, ONE, 0, t[1], 0, 0, ONE, t[2], 0, 0, 0, ONE)
def mat_scale(s):
"""缩放:对角元为 Q16.16 缩放因子。"""
return (s[0], 0, 0, 0, 0, s[1], 0, 0, 0, 0, s[2], 0, 0, 0, 0, ONE)
def mat_rot_y(c, s):
"""绕 Y 轴旋转:c=cos, s=sin(Q16.16);c=0,s=ONE 即 90 度,逐项精确。"""
return (c, 0, s, 0, 0, ONE, 0, 0, -s, 0, c, 0, 0, 0, 0, ONE)
def max_intermediate(a, b):
"""矩阵乘 acc=sum_k a[r,k]*b[k,c] 的最大 |单项乘积| 与最大 |累加|。"""
mp = ma = 0
for r in range(4):
for c in range(4):
acc = 0
for k in range(4):
p = a[r * 4 + k] * b[k * 4 + c]
mp = max(mp, abs(p))
acc += p
ma = max(ma, abs(acc))
return mp, ma
def main():
S = mat_scale((2 * ONE, 2 * ONE, 2 * ONE))
R = mat_rot_y(0, ONE) # 90 度
T = mat_translate((ONE, -ONE, 2 * ONE))
M = cglib.mat_mul(T, cglib.mat_mul(R, S)) # 先缩放 -> 再旋转 -> 再平移
p = (ONE, 2 * ONE, 3 * ONE)
q = cglib.mat_vec4(M, (p[0], p[1], p[2], ONE))
step = cglib.mat_vec4(T, cglib.mat_vec4(R, cglib.mat_vec4(S, (p[0], p[1], p[2], ONE))))
q_tr = cglib.mat_vec4(cglib.mat_mul(T, R), (p[0], p[1], p[2], ONE))
q_rt = cglib.mat_vec4(cglib.mat_mul(R, T), (p[0], p[1], p[2], ONE))
spec = cglib.load_spec()
cam = spec["camera"]
aspect = cglib.fdiv(cglib.f_from_int(spec["width"]), cglib.f_from_int(spec["height"]))
view = cglib.mat_look_at(tuple(cam["eyeQ16"]), tuple(cam["targetQ16"]), tuple(cam["upQ16"]))
proj = cglib.mat_perspective(cam["tanHalfFovYQ16"], aspect, cam["nearQ16"], cam["farQ16"])
mp, ma = max_intermediate(proj, view)
def f3(v):
return "(" + ", ".join(f"{x / ONE:+.3f}" for x in v[:3]) + ")"
print(f"point (world) : {f3(p)}")
print(f"M·p : {f3(q)} raw x={q[0]}")
print(f"step-by-step : {f3(step)} raw x={step[0]}")
print(f"exact match : {list(q) == list(step)}")
print(f"T·R·p : {f3(q_tr)}")
print(f"R·T·p : {f3(q_rt)}")
print(f"commutative : {list(q_tr) == list(q_rt)}")
print(f"max |a·b| proj·view : {mp} ≈ 2^{math.log2(mp):.2f}")
print(f"max |acc| proj·view : {ma} ≈ 2^{math.log2(ma):.2f}")
print(f"int32 signed max : {2**31 - 1} = 2^{math.log2(2**31 - 1):.2f}")
return 0
if __name__ == "__main__":
sys.exit(main())
实跑输出:
point (world) : (+1.000, +2.000, +3.000)
M·p : (+7.000, +3.000, +0.000) raw x=458752
step-by-step : (+7.000, +3.000, +0.000) raw x=458752
exact match : True
T·R·p : (+4.000, +1.000, +1.000)
R·T·p : (+5.000, +1.000, -2.000)
commutative : False
max |a·b| proj·view : 107580544800 ≈ 2^36.65
max |acc| proj·view : 105912729600 ≈ 2^36.62
int32 signed max : 2147483647 = 2^31.00
三点读数。第一,M·p 与 step-by-step 逐项相等(exact match: True)——本算例里缩放因子是 2、旋转是 90 度、平移是整数,所有中间乘积都是 2 的幂或整数,>> 16 不产生任何截断,于是复合与分步在整数下完全一致。第二,T·R·p = (4, 1, 1) 而 R·T·p = (5, 1, -2),commutative: False 直接给出了「顺序敏感」的数值证据。第三,量级一栏留到下一节展开。
量级:Q16.16 矩阵乘为什么要 BigInt
矩阵乘的每一项是 acc = Σ_k a[r,k]·b[k,c],先累加再右移。两个 Q16.16 数相乘,乘积的量级是两者量级之和:元素若到 2^17,单项乘积就到 2^34,累加 4 项更上一层。实测 mat_mul(proj, view) 的最大单项乘积是 107580544800 ≈ 2^36.65,最大累加是 105912729600 ≈ 2^36.62——比 32 位有符号整数的上限 2^31 - 1 = 2147483647 高出约 50 倍。
这正是 JavaScript 侧必须用 BigInt 的原因。Number 的 >> 运算符在规范里先做 ToInt32,只对低 32 位有效:一个 2^36.65 量级的中间量喂进去,高位被静默丢弃,结果错了却不报错。Python 的 int 是任意精度,BigInt 也是,二者在「先乘后移」上语义一致,才谈得上跨语言逐字节相同。逐像素层则相反——屏幕坐标 Q6、边函数、插值分子的量级已证 ≤ 2^49.6 < 2^53,落在 Number 的精确整数区内,所以逐像素路径用 Number 就够,不必背 BigInt 的开销。这条「顶点层 BigInt、像素层 Number」的分界,是量级算出来的,不是风格偏好。
P 层:跨语言逐项一致
TypeScript 侧用 matMul / matVec4 复现同一批整数,并把 S 侧产出的金标准写死在断言里。最硬的一条是量级演示:同一个中间量,BigInt 右移得 1641548,而 Number 右移因 32 位截断只得 3148——两者必须分叉,这才证明「用 BigInt」不是多余。
// file: lantern/core/test/cg_raster_03_1_model_matrix.test.ts
// tier: P
// verify: node --test 'core/test/cg_raster_03_1_model_matrix.test.ts'
// 3.1 模型矩阵的 P 层:与 S 层 sim/cg_raster_03_1_model_matrix.py 逐项一致。
//
// 金标准整数由 S 侧独立产出:
// /opt/homebrew/bin/python3 sim/cg_raster_03_1_model_matrix.py
// 两侧是两份独立代码(Python int vs BigInt),矩阵乘每项 `>> 16` 后必须相同。
import assert from "node:assert/strict";
import { test } from "node:test";
import { ONE } from "../src/lantern/fixed.ts";
import { matIdentity, matMul, matVec4, type Mat4 } from "../src/lantern/mat.ts";
// ---- 由 S 侧产出的金标准(Q16.16 整数) ----
const GOLDEN_M_TRS: Mat4 = [
0n, 0n, 131072n, 65536n,
0n, 131072n, 0n, -65536n,
-131072n, 0n, 0n, 131072n,
0n, 0n, 0n, 65536n,
];
const GOLDEN_M_DOT_P = [458752n, 196608n, 0n, 65536n];
const GOLDEN_TR_DOT_P = [262144n, 65536n, 65536n, 65536n];
const GOLDEN_RT_DOT_P = [327680n, 65536n, -131072n, 65536n];
// S 侧实测:mat_mul(proj, view) 的最大单项乘积
const GOLDEN_MAX_PRODUCT = 107580544800n;
const S: Mat4 = [2n * ONE, 0n, 0n, 0n, 0n, 2n * ONE, 0n, 0n, 0n, 0n, 2n * ONE, 0n, 0n, 0n, 0n, ONE];
const R: Mat4 = [0n, 0n, ONE, 0n, 0n, ONE, 0n, 0n, -ONE, 0n, 0n, 0n, 0n, 0n, 0n, ONE];
const T: Mat4 = [ONE, 0n, 0n, ONE, 0n, ONE, 0n, -ONE, 0n, 0n, ONE, 2n * ONE, 0n, 0n, 0n, ONE];
const P: [bigint, bigint, bigint, bigint] = [ONE, 2n * ONE, 3n * ONE, ONE];
test("复合矩阵 M = T·R·S 与 S 层逐项相同", () => {
const M = matMul(T, matMul(R, S));
assert.deepEqual(M, GOLDEN_M_TRS);
});
test("M·p 与 S 层相同,且与逐点手算(S->R->T)逐项一致", () => {
const M = matMul(T, matMul(R, S));
const q = matVec4(M, P);
assert.deepEqual(q, GOLDEN_M_DOT_P);
let step = matVec4(S, P);
step = matVec4(R, step);
step = matVec4(T, step);
assert.deepEqual(step, q, "定点矩阵链在此取值下逐项精确,复合与分步必须相等");
});
test("复合不可交换:T·R·p 与 R·T·p 不同", () => {
const qTR = matVec4(matMul(T, R), P);
const qRT = matVec4(matMul(R, T), P);
assert.deepEqual(qTR, GOLDEN_TR_DOT_P);
assert.deepEqual(qRT, GOLDEN_RT_DOT_P);
assert.notDeepEqual(qTR, qRT);
});
test("量级:矩阵乘中间量超 32 位,Number 的 >> 会静默出错", () => {
// Number 的 `>>` 先做 ToInt32,只对低 32 位有效;BigInt 是无限精度。
const big = GOLDEN_MAX_PRODUCT >> 16n;
const num = Number(GOLDEN_MAX_PRODUCT) >> 16; // 32 位截断
assert.notEqual(BigInt(num), big, "32 位 >> 与 BigInt >> 结果必须分叉");
assert.equal(big, 1641548n);
assert.equal(num, 3148);
});
test("单位阵恒等:matMul(I, M) == M", () => {
const M = matMul(T, matMul(R, S));
assert.deepEqual(matMul(matIdentity(), M), M);
});
$ node --test 'core/test/cg_raster_03_1_model_matrix.test.ts'
✔ 复合矩阵 M = T·R·S 与 S 层逐项相同
✔ M·p 与 S 层相同,且与逐点手算(S->R->T)逐项一致
✔ 复合不可交换:T·R·p 与 R·T·p 不同
✔ 量级:矩阵乘中间量超 32 位,Number 的 >> 会静默出错
✔ 单位阵恒等:matMul(I, M) == M
ℹ tests 5 pass 5 fail 0
小结
模型矩阵把局部坐标搬到世界坐标,由平移、旋转、缩放三种 4x4 齐次矩阵复合而成。复合顺序不可交换,T·R 与 R·T 作用到同一点给出不同结果;矩阵作用在右侧,所以 T·R·S 的求值顺序是「先缩放、再旋转、最后平移」。全整数实现下,只要各项乘积不触发截断,复合矩阵与逐点分步就逐项精确一致。量级上,矩阵乘的中间量可达 2^36.65,远超 32 位,故 Python 用 int、TypeScript 必须用 BigInt;而逐像素层量级 ≤ 2^49.6 < 2^53,用 Number 即可。下一节把相机本身也写成矩阵,构造视图矩阵。
阅读导航:上一节:2.3 齐次坐标与投影 · 下一节:3.2 视图矩阵 。
继续阅读
探索更多技术文章
浏览归档,发现更多关于系统设计、工具链和工程实践的内容。