3.1 模型矩阵

从平移、旋转、缩放的 4x4 齐次形式出发,解释复合变换为何不可交换,给出「先缩放、再旋转、再平移」的 Q16.16 定点实现并断言它与逐点手算逐项一致,再用实测数据说明矩阵乘的中间量为何超过 32 位、JavaScript 侧必须改用 BigInt。

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 视图矩阵 。

继续阅读

探索更多技术文章

浏览归档,发现更多关于系统设计、工具链和工程实践的内容。

全部文章 返回首页

「cg」更多文章

  1. 14.3 中间缓冲的可视化
  2. 14.2 覆盖率与线框
  3. 14.1 深度图与法线图