8.3 各向异性过滤
8.2 的 mip 链是各向同性的：每一级在两个纹理方向上都减半。这在足迹接近正方形时没问题，但真实场景里足迹经常是细长条——掠射角的地面、斜看的墙面、伸向远方的道路，屏幕像素在纹理空间里跨的 $u$ 和跨的 $v$ 能差一个数量级。本节处理这种情形。边界照旧：纹理与显存 那篇管纹理的存储与驻留，本节管的是过滤核的形状——足迹怎么算、金字塔怎么改成两维、采样怎么近似。
各向同性过滤的病：细长足迹
8.2 选级用的足迹是 $\max(|\partial u/\partial x|, |\partial v/\partial y|) \cdot n$。这个 $\max$ 是问题的根源。
设想一个像素在 $u$ 方向跨了 16 个 texel、在 $v$ 方向只跨了 1 个。取 $\max$ 得到 16，于是 LOD 选到第 4 级——$16 \times 16$ 的块被平均成 1 个 texel。可是 $v$ 方向本来只需要跨 1 个 texel，这一平均把 $v$ 方向也糊掉了 16 倍。画面上的表现是：远处的路面不但没有闪烁，反而糊成一片，细节全没了。
反过来取 $\min$ 呢？选到第 0 级，$u$ 方向那 16 个 texel 只采了 2 个（双线性），剩下 14 个全丢——闪烁回来了。
两个方向需要不同的缩放，而 mip 链只提供一个数字。 这就是各向异性过滤要解决的问题。
从 UV 梯度到椭圆足迹
屏幕像素是一个正方形，它在 UV 空间的像是一个平行四边形，由雅可比矩阵给出：
$$J = \begin{bmatrix} \partial u/\partial x & \partial u/\partial y \ \partial v/\partial x & \partial v/\partial y \end{bmatrix}$$
四行梯度都是光栅化时顺手算出来的（透视校正插值本来就要它们的倒数），不需要额外采样。过滤核取这个平行四边形的内切椭圆——只有椭圆的半轴才给出「一个方向缩多少、另一个方向缩多少」。半轴平方就是 $J J^{\mathsf{T}}$ 的两个特征值：
$$a = \Big(\tfrac{\partial u}{\partial x}\Big)^2 + \Big(\tfrac{\partial u}{\partial y}\Big)^2, \quad
c = \Big(\tfrac{\partial v}{\partial x}\Big)^2 + \Big(\tfrac{\partial v}{\partial y}\Big)^2, \quad
b = \tfrac{\partial u}{\partial x}\tfrac{\partial v}{\partial x} + \tfrac{\partial u}{\partial y}\tfrac{\partial v}{\partial y}$$
$$\lambda_{\max,\min} = \frac{a + c}{2} \pm \sqrt{\Big(\frac{a - c}{2}\Big)^2 + b^2}, \qquad
\text{major} = \sqrt{\lambda_{\max}}, \quad \text{minor} = \sqrt{\lambda_{\min}}$$
长轴方向取 $\lambda_{\max}$ 的特征向量，对 $2 \times 2$ 对称阵可以直接写出：$(\lambda_{\max} - c,; b)$，归一化即可。$\text{major} / \text{minor}$ 就是各向异性比。
自检里的两个例子把这件事说清楚了：梯度 $(1, 0)$ 配 $(0, 1/8)$（$u$ 方向每像素走 1 个 uv 单位、$v$ 方向走 $1/8$）给出 major 65536、minor 8192、比 8，长轴方向 $(65536, 0)$；两个方向都走 1 时 major = minor = 65536、比 1。比 1 的像素不需要各向异性——而比大于 1 的像素有多少，才是这一节值不值得做的判据。
注意这里有个纯定点的坑：$a$、$c$、$b$ 都是 Q16.16 乘法的结果，而梯度本身可能很小（掠射角下 $\partial u/\partial x$ 只有几十）。fmul(x, x) = (x * x) >> 16 在 $x < 256$ 时直接把平方截成 0。这不是 bug，是定点分辨率的必然代价；本节场景的梯度足够大（最远处约 $5.5 \times 10^3$），所以没有踩到，但把它写在这里，是因为换一个更远的场景就会踩。
同一段公式搬到 TypeScript 侧还有第二个坑：$\frac{a - c}{2}$ 和 $\frac{a + c}{2}$ 都不能写成 >> 1。$a$、$c$ 是 Q16.16 的乘积，量级在 $2^{32}$ 上下，而 JS 的移位只对 32 位有效，一移位就变成负数——必须显式写成除法。下面是 cg_shade_08_aniso.ts 里 ellipse 的前半段，逐字摘录：
// snippet: lantern/core/src/lantern/cg_shade_08_aniso.ts（`ellipse` 前半段，逐字）
const a = fmul(gx[0], gx[0]) + fmul(gx[1], gx[1]);
const c = fmul(gy[0], gy[0]) + fmul(gy[1], gy[1]);
const b = fmul(gx[0], gy[0]) + fmul(gx[1], gy[1]);
const d0 = Math.floor((a - c) / 2); // 不用 `>> 1`：JS 的移位只对 32 位有效
const d1 = b;
const h = fsqrt(fmul(d0, d0) + fmul(d1, d1));
const half = Math.floor((a + c) / 2);
const lmax = half + h;
const lmin = half - h;
const major = fsqrt(Math.max(lmax, 0));
const minor = fsqrt(Math.max(lmin, 0));
Python 侧因为 int 是无限精度，同一段直接写 (a - c) >> 1 就对。这正是「先跑 S 层内核、再移植到 P 层」这套顺序能提前抓住的差异：cg_shade_08_3_aniso.test.ts 里那条「定点原语与 Python 一致」的用例就是钉这个的，该文件 6 条用例全绿。
Ripmap：两个方向各自索引
既然两个方向需要不同的缩放，最直接的办法是把金字塔做成两维：ripmap，$L(i, j)$ 表示 $u$ 方向缩 $2^i$、$v$ 方向缩 $2^j$。构造分两步——先沿 $u$ 折行造出 $L(i, 0)$，再对每个 $i$ 沿 $v$ 折列造出 $L(i, j)$：
$$L(i, 0) = \big(L(i-1, 0)[:, 0::2] + L(i-1, 0)[:, 1::2]\big) \gg 1, \qquad
L(i, j) = \big(L(i, j-1)[0::2, :] + L(i, j-1)[1::2, :]\big) \gg 1$$
$32 \times 32$ 的纹理造出 $6 \times 6 = \mathbf{36}$ 级，每个像素一次查表就能拿到「两个方向各自合适」的那一级，走样量是全场最低的 11.00。
代价是显存。总 texel 数是
$$\sum_{i=0}^{\text{lv}} \sum_{j=0}^{\text{lv}} \lfloor n/2^i \rfloor \lfloor n/2^j \rfloor = \Big(\sum_{i=0}^{\text{lv}} \lfloor n/2^i \rfloor\Big)^2$$
$n = 32$ 时内层和是 $32+16+8+4+2+1 = 63$，平方得 3969 texel，相对原始的 $1024$ 是 ×3.88（3875‰）。渐近界是 ×4：$n \to \infty$ 时内层和趋于 $2n$，比值趋于 $(2n)^2 / n^2 = 4$，本节量到的 3.88 正在往它靠。作为对照，各向同性的 mip 链只占 $1 + 1/4 + 1/16 + \cdots = 4/3$，$32 \times 32$ 上是 $1365$ texel、×1.33——ripmap 比它贵大约三倍。很多资料里那句「mip 链只多占三分之一显存」说的正是各向同性链，不是 ripmap；引用之前先算一遍自己的尺寸。
N 采样：真机的近似
存整张 ripmap 太贵，真机用的是各向异性 N 采样：层级按短轴取（用普通的各向同性链就够），再沿长轴取 $N$ 个样本平均。样本沿长轴均匀铺开：
$$t_k = \frac{(2k + 1 - N) \cdot \text{major}}{2N}, \qquad k = 0, 1, \dots, N-1$$
即 $t$ 关于中心对称、间隔 $\text{major}/N$，覆盖 $[-\text{major}/2, +\text{major}/2)$。样本数不必等于 $N$——如果足迹本身就是圆的（major ≈ minor），一个样本就够，所以实际取 $\min\big(N, \lceil \text{major} / \text{minor} \rceil\big)$。这一点很重要：只有各向异性比超过 1 的像素才付出多采样的代价。
内核（S 层）
本节的内核落在 sim/cg_shade_08_3_kernel.py：噪声、各向同性链、ripmap、非方阵双线性、$\lfloor \log_2 \rfloor$、以及把梯度变成椭圆的那个函数。驱动脚本负责场景、参考解与曲线。
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_08_3_kernel.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_08_3_kernel.py
"""8.3 的内核：UV 梯度 -> 椭圆足迹 -> Ripmap -> 各向异性 N 采样。
被 sim/cg_shade_08_3_aniso.py 调用；单独成文件是为了让正文里的 `file:` 块就是
本节真正在跑的那段代码。
屏幕像素正方形在 UV 空间的像是平行四边形，过滤核取它的**内切椭圆**——只有椭圆的
半轴才给出「一个方向缩多少、另一个方向缩多少」。设雅可比
J = [[du/dx, du/dy], [dv/dx, dv/dy]] （单位：texel / 像素）
则半轴平方是 `J Jᵀ` 的特征值 `λ = (a+c)/2 ± sqrt(((a-c)/2)² + b²)`，其中
`a = du/dx² + du/dy²`、`c = dv/dx² + dv/dy²`、`b = du/dx·dv/dx + du/dy·dv/dy`。
长轴方向取 lmax 的特征向量 `(lmax - c, b)`，半轴比就是各向异性比。
**Ripmap** = 按 (i, j) 两维索引的 mip 金字塔，L(i,j) 在 u 方向缩 2^i、v 方向
缩 2^j。单次查表即可拿到「两个方向各自合适」的那一级。本文件不含浮点。
"""
import math
import os
import sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ONE = cglib.ONE
def make_noise(n, seed=0x0C7A):
"""程序化高频纹理：每 2x2 texel 共享一个整数哈希值（`cglib._hash_noise`）。"""
tex = np.zeros((n, n, 3), dtype=np.int64)
for ty in range(n):
for tx in range(n):
for ch in range(3):
tex[ty, tx, ch] = cglib._hash_noise(tx >> 1, ty >> 1, seed + ch) >> 8
return tex
def build_mips(tex):
"""各向同性 mip 链：两个方向同时减半，`(a+b+c+d) >> 2`。"""
chain, cur = [tex], tex
while cur.shape[0] > 1:
cur = (cur[0::2, 0::2] + cur[0::2, 1::2] + cur[1::2, 0::2] + cur[1::2, 1::2]) >> 2
chain.append(cur)
return chain
def build_ripmap(tex):
"""Ripmap：L(i,j) 在 u 方向缩 2^i、v 方向缩 2^j。先沿 u 折行，再沿 v 折列。"""
lv = max(0, int(math.log2(tex.shape[0])))
rip = [[None] * (lv + 1) for _ in range(lv + 1)]
rip[0][0] = tex
for i in range(1, lv + 1):
rip[i][0] = (rip[i - 1][0][:, 0::2] + rip[i - 1][0][:, 1::2]) >> 1
for i in range(lv + 1):
for j in range(1, lv + 1):
rip[i][j] = (rip[i][j - 1][0::2] + rip[i][j - 1][1::2]) >> 1
return rip
def sample_bilinear(tex, u, v):
"""双线性：texel 中心在 (i+0.5)/n；ripmap 的级非方阵，故行数列数分开取。"""
ny, nx = tex.shape[0], tex.shape[1]
fu, fv = u * nx - (ONE >> 1), v * ny - (ONE >> 1)
iu, iv = fu >> 16, fv >> 16
gu, gv = fu - (iu << 16), fv - (iv << 16)
a, b, c, d = ONE - gu, gu, ONE - gv, gv
x0, x1 = iu % nx, (iu + 1) % nx
y0, y1 = iv % ny, (iv + 1) % ny
return [(int(tex[y0, x0, ch]) * a * c + int(tex[y0, x1, ch]) * b * c
+ int(tex[y1, x0, ch]) * a * d + int(tex[y1, x1, ch]) * b * d) >> 32
for ch in range(3)]
def ilog2(x):
"""floor(log2(x))：整数位长求法，不用浮点 log。"""
return int(x).bit_length() - 1
def ellipse(gx, gy):
"""由 (du/dx, dv/dx, du/dy, dv/dy) 的 Q16.16 梯度求椭圆半轴与长轴方向（uv 单位）。"""
a = cglib.fmul(gx[0], gx[0]) + cglib.fmul(gx[1], gx[1])
c = cglib.fmul(gy[0], gy[0]) + cglib.fmul(gy[1], gy[1])
b = cglib.fmul(gx[0], gy[0]) + cglib.fmul(gx[1], gy[1])
d0, d1 = (a - c) >> 1, b
h = cglib.fsqrt(cglib.fmul(d0, d0) + cglib.fmul(d1, d1))
lmax, lmin = ((a + c) >> 1) + h, ((a + c) >> 1) - h
major, minor = cglib.fsqrt(max(lmax, 0)), cglib.fsqrt(max(lmin, 0))
if b == 0: # 轴对齐：长轴就是 u 或 v 方向
ex, ey = (ONE, 0) if a >= c else (0, ONE)
else: # J Jᵀ 的 lmax 特征向量 = (lmax - c, b)
ex, ey = lmax - c, b
n = cglib.fsqrt(cglib.fmul(ex, ex) + cglib.fmul(ey, ey))
ex, ey = (cglib.fdiv(ex, n), cglib.fdiv(ey, n)) if n > 0 else (ONE, 0)
return major, minor, ex, ey
def level_of(footprint):
"""texel 单位的轴足迹 -> mip 级号：floor(log2(footprint)) - 16，钳到 [0, lv]。"""
return max(0, ilog2(max(1, footprint)) - 16)
if __name__ == "__main__":
t = make_noise(16)
rip = build_ripmap(t)
lv = len(rip) - 1
total = sum(int(rip[i][j].shape[0]) * int(rip[i][j].shape[1])
for i in range(lv + 1) for j in range(lv + 1))
print("16x16 ripmap :", lv + 1, "级", 16 * 16, "->", total, "texel")
print("各向同性链尺寸 :", [int(l.shape[0]) for l in build_mips(t)])
print("L(0,3) 尺寸 :", rip[0][3].shape, " L(3,0) 尺寸:", rip[3][0].shape)
ma, mi, ex, ey = ellipse((ONE, 0), (0, ONE >> 3))
print("u 走 1 / v 走 1/8 : major", ma, "minor", mi, "比", ma // mi, "长轴", (ex, ey))
ma, mi, ex, ey = ellipse((ONE, 0), (0, ONE))
print("两个方向都走 1 : major", ma, "minor", mi, "比", ma // mi)
print("level_of(1 texel) :", level_of(ONE), " (8 texel):", level_of(8 * ONE))
自检输出（/opt/homebrew/bin/python3 sim/cg_shade_08_3_kernel.py）：
16x16 ripmap : 5 级 256 -> 961 texel
各向同性链尺寸 : [16, 8, 4, 2, 1]
L(0,3) 尺寸 : (2, 16, 3) L(3,0) 尺寸: (16, 2, 3)
u 走 1 / v 走 1/8 : major 65536 minor 8192 比 8 长轴 (65536, 0)
两个方向都走 1 : major 65536 minor 65536 比 1
level_of(1 texel) : 0 (8 texel): 3
$16 \times 16$ 那张图的 ripmap 是 $31^2 = 961$ texel，正好是内层和的平方——同一个公式换一个 $n$ 而已。注意 L(0,3) 的形状是 (2, 16, 3)（2 行、16 列）而 L(3,0) 是 (16, 2, 3)：ripmap 的级不是方阵，所以采样时行数列数必须分开取，不能只用一个 $n$。
走样曲线与负例
场景换成一块掠射角地面（$96 \times 72$，四边形从 $z = 14$ 伸到 $z = -28$），覆盖 984 像素。各向异性比的分布是 [377, 222, 273, 111, 0, 1]（按 $2^k$ 分档）：比值小于 2 的只有 377 个，其余 607 个（62%）都超过 2，还有一个超过 32。这块场景确实需要各向异性。
走样量同样对 $8 \times 8$ 超采样参考解算 MAE：
| 模式 | MAE ↓ | 耗时 | 相对上一档 |
|---|---|---|---|
| $N = 1$ | 25.9915 | 4.0 ms | — |
| $N = 2$ | 24.0766 | 5.7 ms | $-7.4%$ |
| $N = 4$ | 20.5301 | 8.8 ms | $-14.7%$ |
| $N = 8$ | 19.7619 | 11.1 ms | $-3.7%$ |
| ripmap | 11.0030 | 3.3 ms | — |
rip-major（负例） | 14.3496 | 3.4 ms | $+30.4%$ |
rip-minor（负例） | 25.9959 | 4.5 ms | $+136.3%$ |
| 参考解（$8 \times 8$） | — | 272.3 ms | — |
三条结论。
一、样本数买到的收益是递减的。 $N$ 从 1 到 2 改善 7.4%，2 到 4 改善 14.7%，4 到 8 只剩 3.7%——而耗时几乎线性（4.0 → 5.7 → 8.8 → 11.1 ms）。把 $N$ 从 4 提到 8，多花约三成的时间换 3.7% 的质量，这笔账在实时渲染里通常是划不来的。默认值取 4 或 8 是有道理的，再往上就不值了。
二、ripmap 又便宜又好，但它要显存。 单次查表 3.3 ms、MAE 11.00，比 $N = 8$ 的 11.1 ms / 19.76 好得多——因为 N 采样只是用有限个样本去近似椭圆，而 ripmap 是精确的两个方向各自缩放。代价是 ×3.88 的显存，在移动端基本不可接受。所以真实硬件走 N 采样，ripmap 主要出现在离线渲染和对照实验里。
三、两个索引必须独立，绑到对角线上立刻退化。 这是本节留的两条负例：把 ripmap 的两个层级索引强行绑到对角线，rip-minor（两轴都取较细的那级）的 MAE 是 25.9959——和 $N = 1$（25.9915）只差 0.0044，也就是说它和「完全不看各向异性、直接用各向同性链」几乎完全等价；rip-major（都取较粗的那级）是 14.3496，比正确的 ripmap 差 30.4%，属于过度模糊。两条都远超 1.0 的容差判据。
跨宿主断言（P 层）
TypeScript 侧的 core/src/lantern/cg_shade_08_aniso.ts 与内核逐字对齐。它比 8.2 多一个原语要小心：$\sqrt{\cdot}$ 用的是整数平方根 isqrt(x << 16)，而 $x \ll 16$ 在掠射角下会逼近 $2^{53}$，double 的 Math.sqrt 在那里不再可靠——所以 TS 侧用 BigInt 牛顿法，保证逐位一致。
✔ 定点原语：isqrt / fsqrt 与 Python 的 math.isqrt 一致
✔ 【核心】ripmap 逐级 sha256 与 Python 侧逐字节相同
✔ 【核心】覆盖数与各向异性比直方图一致
✔ 【核心】N=1/2/4/8 与三种 ripmap 规则的像素缓冲 sha256 逐字节相同
✔ 8x8 超采样参考解逐字节相同，且 N 曲线单调下降
✔ 负例：ripmap 把两个层级索引绑到对角线必须被 1.0 的容差判据抓住
ℹ tests 6 ℹ pass 6 ℹ fail 0
36 个 ripmap 级、$N$ 的四种取值、三种 ripmap 规则、加一张参考解——共 44 张像素缓冲的 sha256 全部与 Python 侧逐字节相同。其中「N 曲线单调下降」也是一条断言：如果哪次改动让 $N = 8$ 比 $N = 4$ 更差，说明采样方向或层级选择写反了。
小结
各向异性过滤把 8.2 里那个 $\max$ 拆开了：足迹从「一个数字」变成「一个椭圆」，金字塔从一维变成两维，选级从「一个索引」变成「两个索引」。走样量从 25.99 压到 11.00（ripmap）或 19.76（$N = 8$ 采样），代价分别是 ×3.88 的显存和 ×2.8 的耗时。到这里，纹理这条线——坐标、寻址、过滤、金字塔、各向异性——都落在整数域、都能跨宿主逐位断言了。下一章换一个方向：不再从纹理里查颜色，而是用坐标本身算出颜色，也就是过程化噪声。
阅读导航：上一节：8.2 过滤与 mipmap 的构造 · 下一节：9.1 值噪声与 Perlin 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。