9.2 分形噪声与 fBm
站内 /graphics-terrain-large-world/
那篇关心的是把分形地形切成块、按距离换 LOD、再流式塞进显存；它默认「高度场已经有了」。本节补上那一步：高度场本身怎么由整数噪声叠出来。单层噪声太规整——每个格点一个值，一眼就能看出 4 单位的网格。分形噪声（fBm，fractional Brownian motion）把不同尺度的噪声叠在一起，用三个旋钮控制「粗糙」的程度。三个旋钮全在整数域：octaves 是层数，lacunarity 是每层格点间距的倍率（这里取 2，由 shift - o 直接实现），gain 是每层幅度的衰减（这里 1/2，fmul(amp, gain)）。没有浮点乘，也没有除法——倍频是 2 的幂，直接作用在移位量上。
整数 fBm 的三行核心
{
"name": "lantern-fbm-scan",
"grid": 128,
"stepQ16": 8192,
"latticeShift": 18,
"seed": 24285,
"gainQ16": 32768,
"octaves": [1, 2, 4, 8],
"spectrumBins": 48,
"warpOctaves": 4,
"warpStrengthQ16": 131072,
"probes": [[0, 0], [1, 0], [2, 3], [5, 7], [11, 13], [31, 17], [64, 64], [100, 50],
[127, 127], [40, 80]]
}
叠加循环只有三行，值得逐行看：
# snippet: 节选自 lantern/sim/cg_shade_09_noise.py 的 fbm() / warp()
def fbm(fx, fz, seed, octaves=4, shift=18, gain=ONE >> 1):
total, amp, norm = 0, ONE, 0
for o in range(octaves):
n = value_noise(fx, fz, seed + o, shift - o) # 格点间距每层减半
total += cglib.fmul(n, amp) # 幅度每层乘 gain
norm += amp
amp = cglib.fmul(amp, gain)
return cglib.fdiv(total, norm) # 分子非负，fdiv == floor
def warp(fx, fz, seed, strength, octaves=4, shift=18):
ox = fbm(fx, fz, seed + 101, octaves, shift) - (ONE >> 1) # 居中，有正有负
oz = fbm(fx, fz, seed + 202, octaves, shift) - (ONE >> 1)
return fx + cglib.fmul(ox, strength), fz + cglib.fmul(oz, strength)
第一行 value_noise(fx, fz, seed + o, shift - o) 同时做了两件事：shift - o 让格点间距每层减半（lacunarity = 2），seed + o 让每层用不同的哈希，否则各层会在格点上对齐、叠加出可见的网格。第二行 total += fmul(n, amp) 是幅度加权。第三行 fdiv(total, norm) 按幅度和归一化回 [0, ONE)——这一步不能省，否则层数一多，和就会溢出 ONE，颜色和地形高度都得再夹一次。
fdiv 在这里是安全的：值噪声恒非负，total 与 norm 都非负，所以「向零截断」与「向下取整」是同一个结果。这也是 BRIEF 铁律「降精度用 >>、除法只用在已知非负的分子上」的一个具体落点。
域扭曲（domain warping）是同一套东西的二次消费：先用两个低频 fBm 算出一个偏移量，再拿偏移后的坐标去采样。ox、oz 减掉 ONE/2 是为了居中——偏移量必须有正有负，否则整张图会被整体推走而不是被「拧」。同样没有除法，fmul(ox, strength) 直接缩放。
扫描 octaves = 1 / 2 / 4 / 8
用 0.125 世界单位步长（每基准格 32 点）采样 128×128，让最细的第 8 层（格点间距 0.03125）也被充分采样，否则量到的是混叠而不是噪声：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_shade_09_2_fbm.py
网格 128x128 步长 0.125 基准格点 4.0 gain 0.500
octaves mean_q16 std min max roughness rough/std centroid ms
1 25624 13653.3 169 57539 553.6 0.0405 0.0067 41.5
2 27516 10845.1 4438 58662 552.3 0.0509 0.0073 89.4
4 28563 8995.2 5510 53406 642.5 0.0714 0.0074 155.2
8 28823 8481.3 6951 53028 778.6 0.0918 0.0074 316.3
域扭曲（o=4, 强度 2.0）：与 o=4 基准 |Δ| 均值 1571 归一化粗糙度 0.0732（未扭曲 0.0714） 耗时 480.8 ms
方差随层数下降：13653 → 10845 → 8995 → 8481。 这违反「细节越多越花」的直觉。原因是归一化：norm 是幅度和，octaves=1 时是 ONE，octaves=8 时约 2×ONE。低频那一层（幅度最大、方差最大）在总和中被稀释了，所以总方差反而降。归一化 fBm 的方差不是层数的单调函数。
归一化粗糙度单调上升：0.0405 → 0.0509 → 0.0714 → 0.0918。 这是本节的判据。它用「相邻像素 |差| 均值」除以「标准差」，把幅度衰减的干扰除掉，只剩「每步起伏相对于总振幅有多大」。8 层比 1 层粗糙 2.27 倍，且严格单调——这才说明加进去的高频层真的落在可分辨的尺度上，而不是被归一化吃掉了。
绝对粗糙度（553.6 → 778.6）反而几乎不单调（2 层 552.3 比 1 层还小）。这正说明单看绝对量会误判：层数增加时低频被稀释、高频被加入，两股力量在绝对量上互相抵消。判据必须归一化，否则量到的是归一化系数而不是噪声的尺度结构。功率谱质心几乎不变（0.0067 → 0.0074），说明 0.125 的采样步长已经接近这张网格能分辨的极限：再加层，能量只能堆在同一个 bin 里。这也提示——质心是网格分辨率的函数，粗糙度比值不是，所以本节把前者当辅助量、后者当主判据。
耗时线性增长：41.5 → 89.4 → 155.2 → 316.3 ms。 每一层都要多跑一遍 value_noise，8 层是 1 层的 7.6 倍。没有捷径：fBm 的代价就是层数乘以单层代价。
域扭曲：把细节拧成条带
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_09_2_fbm.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_09_2_fbm.py
"""9.2 分形噪声与 fBm：倍频叠加的整数实现、粗糙度扫描与域扭曲。
fBm = 逐层叠加同一个 `value_noise`：每层格点间距减半（lacunarity = 2，由 `shift - o`
直接实现，不需要浮点乘），幅度乘以 `gain`（这里 1/2），最后按幅度和归一化。三个旋钮
`octaves` / `lacunarity` / `gain` 全在整数域。
扫描 octaves = 1 / 2 / 4 / 8，用 **0.125 世界单位步长**（每基准格 32 点）量三件事：
* 相邻像素 |差| 均值（绝对粗糙度）；
* 归一化粗糙度 = |Δ| 均值 / 标准差（去掉幅度衰减的干扰，看「每步起伏」）；
* 功率谱质心（能量集中在多高的空间频率）。
再做域扭曲（domain warping）：用两个低频 fBm 偏移采样坐标，细节被拧成条带。
产物：out/cg_shade_09_2_fbm.{csv,json}、out/cg_shade_09_2_fbm_o{1,2,4,8}.png、
out/cg_shade_09_2_warp.png。
"""
import json
import os
import sys
import time
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cg_shade_09_noise as nz # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
SPEC = os.path.join(ROOT, "config", "cg_shade_09_2_fbm.json")
def sample_grid(n, step, fn):
a = np.zeros((n, n), dtype=np.int64)
for y in range(n):
fz = y * step
for x in range(n):
a[y, x] = fn(x * step, fz)
return a
def roughness(field):
"""相邻像素 |差| 均值（Q16.16）；再除以标准差得到无量纲的「每步起伏」。"""
d = np.concatenate([np.abs(np.diff(field, axis=1)).ravel(),
np.abs(np.diff(field, axis=0)).ravel()])
std = float(field.std())
return float(d.mean()), (float(d.mean()) / std if std > 0 else 0.0)
def centroid(field, bins):
"""功率谱质心（归一化到 [0,1)）：能量越偏高频，值越大。"""
f = field.astype(np.float64)
p = np.abs(np.fft.fftshift(np.fft.fft2(f - f.mean()))) ** 2
h, w = p.shape
cy, cx = h // 2, w // 2
yy, xx = np.mgrid[0:h, 0:w]
r = np.sqrt((yy - cy) ** 2 + (xx - cx) ** 2)
rmax = min(cy, cx)
idx = np.clip((r / rmax * bins).astype(np.int64), 0, bins - 1)
s = np.array([p[idx == b].mean() if (idx == b).any() else 0.0 for b in range(bins)])
s = s / s.sum()
return float((s * (np.arange(bins) / bins)).sum())
def main():
spec = json.load(open(SPEC, encoding="utf-8"))
n, step, shift = spec["grid"], spec["stepQ16"], spec["latticeShift"]
seed, gain, sbins = spec["seed"], spec["gainQ16"], spec["spectrumBins"]
octs = spec["octaves"]
rows, fields, times = [], {}, []
for o in octs:
t0 = time.perf_counter()
f = sample_grid(n, step, lambda fx, fz, o=o: nz.fbm(fx, fz, seed, o, shift, gain))
ms = (time.perf_counter() - t0) * 1000.0
rough, rel = roughness(f)
fields[o] = f
times.append(ms)
rows.append({"octaves": o, "mean": float(f.mean()), "std": float(f.std()),
"min": int(f.min()), "max": int(f.max()),
"rough": rough, "roughRel": rel, "centroid": centroid(f, sbins),
"ms": round(ms, 3)})
# 域扭曲：坐标被两个低频 fBm 偏移后再取样，细节沿偏移方向拉成条带
wo, wstr = spec["warpOctaves"], spec["warpStrengthQ16"]
t0 = time.perf_counter()
warped = sample_grid(n, step, lambda fx, fz: nz.fbm(*nz.warp(fx, fz, seed, wstr, wo, shift),
seed, wo, shift, gain))
warp_ms = (time.perf_counter() - t0) * 1000.0
base = fields[wo]
warp_delta = float(np.abs(warped - base).mean())
warp_rough, warp_rel = roughness(warped)
probes = []
for p in spec["probes"]:
fx, fz = p[0] * step, p[1] * step
wx, wz = nz.warp(fx, fz, seed, wstr, wo, shift)
probes.append({"x": p[0], "z": p[1], "fx": fx, "fz": fz,
"fbm4": int(nz.fbm(fx, fz, seed, 4, shift, gain)),
"warpX": int(wx), "warpZ": int(wz),
"warpedFbm": int(nz.fbm(wx, wz, seed, wo, shift, gain))})
out = os.path.join(ROOT, "out")
os.makedirs(out, exist_ok=True)
for o in octs:
cglib.save_gray_png(os.path.join(out, f"cg_shade_09_2_fbm_o{o}.png"), fields[o])
cglib.save_gray_png(os.path.join(out, "cg_shade_09_2_warp.png"), warped)
cglib.write_csv(os.path.join(out, "cg_shade_09_2_fbm.csv"),
["octaves", "mean_q16", "std", "min", "max", "roughness", "rough_rel",
"centroid", "ms"],
[[r["octaves"], int(r["mean"]), f"{r['std']:.1f}", r["min"], r["max"],
f"{r['rough']:.1f}", f"{r['roughRel']:.4f}", f"{r['centroid']:.4f}",
r["ms"]] for r in rows])
stats = {
"grid": n, "stepQ16": step, "latticeShift": shift, "seed": seed, "gainQ16": gain,
"octaves": octs, "spectrumBins": sbins, "scan": rows,
"warp": {"octaves": wo, "strengthQ16": wstr, "ms": round(warp_ms, 3),
"deltaFromBase": warp_delta, "roughness": warp_rough, "roughRel": warp_rel,
"mean": float(warped.mean()), "std": float(warped.std())},
"probes": probes,
}
with open(os.path.join(out, "cg_shade_09_2_fbm.json"), "w", encoding="utf-8") as fh:
json.dump(stats, fh, ensure_ascii=False, indent=2)
print(f"网格 {n}x{n} 步长 {step / cglib.ONE:.3f} 基准格点 {2 ** shift / cglib.ONE:.1f} "
f"gain {gain / cglib.ONE:.3f}")
print("octaves mean_q16 std min max roughness rough/std centroid ms")
for r in rows:
print(f"{r['octaves']:>7d} {int(r['mean']):8d} {r['std']:7.1f} {r['min']:5d} "
f"{r['max']:5d} {r['rough']:9.1f} {r['roughRel']:9.4f} {r['centroid']:8.4f} "
f"{r['ms']:6.1f}")
print(f"域扭曲（o={wo}, 强度 {wstr / cglib.ONE:.1f}）：与 o={wo} 基准 |Δ| 均值 {warp_delta:.0f} "
f"归一化粗糙度 {warp_rel:.4f}（未扭曲 {roughness(base)[1]:.4f}） 耗时 {warp_ms:.1f} ms")
print(f"探针 fbm4（前 3） " + " ".join(f"({p['x']},{p['z']})={p['fbm4']}" for p in probes[:3]))
print(f"探针 warp（前 3） " + " ".join(f"({p['x']},{p['z']})->({p['warpX']},{p['warpZ']})"
for p in probes[:3]))
print("png/csv/json: out/cg_shade_09_2_*.{png,csv,json}")
return 0
if __name__ == "__main__":
sys.exit(main())
在 4 层基准上做强度 2.0 的扭曲，实测与未扭曲版本的 |Δ| 均值是 1571（相对标准差约 9000，即约 17% 的像素被明显改动），归一化粗糙度从 0.0714 微升到 0.0732，耗时 480.8 ms（三倍于单次 4 层 fBm——每像素要多算两个偏移用的 fBm）。一个诚实的观察：域扭曲主要改变的是「细节往哪走」，而不是「有多粗糙」。全局粗糙度几乎不动，但纹理的走向被低频场拽成了漩涡与条带。这也是它常被用来做云、大理石、火焰的原因——需要的是方向性，不是更多的频率。
跨语言核对
P 层必须复现 fBm 探针与扭曲后的坐标：
// file: lantern/core/test/cg_shade_09_2_fbm.test.ts
// tier: P
// verify: node --test 'core/test/cg_shade_09_2_fbm.test.ts'
// 9.2 的跨语言断言：TypeScript 必须逐值复现 Python 侧 sim/cg_shade_09_2_fbm.py
// 的 fBm 探针与域扭曲坐标，并用负例证明「归一化」这一步不能省。
import assert from "node:assert/strict";
import { readFileSync } from "node:fs";
import { dirname, join } from "node:path";
import { test } from "node:test";
import { fileURLToPath } from "node:url";
import { ONE, fmul } from "../src/lantern/fixed.ts";
import { fbm, valueNoise, warp } from "../src/lantern/cg_shade_09_noise.ts";
const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const cfg = JSON.parse(
readFileSync(join(root, "out", "cg_shade_09_2_fbm.json"), "utf8"),
) as {
latticeShift: number;
seed: number;
gainQ16: number;
warp: { octaves: number; strengthQ16: number };
probes: { x: number; z: number; fx: number; fz: number; fbm4: number; warpX: number; warpZ: number }[];
};
const SHIFT = BigInt(cfg.latticeShift);
const SEED = BigInt(cfg.seed);
const GAIN = BigInt(cfg.gainQ16);
const WO = cfg.warp.octaves;
const WSTR = BigInt(cfg.warp.strengthQ16);
test("【核心】fBm 探针与域扭曲坐标逐值一致", () => {
for (const p of cfg.probes) {
const fx = BigInt(p.fx);
const fz = BigInt(p.fz);
assert.equal(fbm(fx, fz, SEED, 4, SHIFT, GAIN), BigInt(p.fbm4), `(${p.x},${p.z}) fbm4`);
const [wx, wz] = warp(fx, fz, SEED, WSTR, WO, SHIFT);
assert.equal(wx, BigInt(p.warpX), `(${p.x},${p.z}) warpX`);
assert.equal(wz, BigInt(p.warpZ), `(${p.x},${p.z}) warpZ`);
}
});
test("单层 fBm 退化为值噪声（归一化系数恰为 ONE）", () => {
for (let k = 0; k < 50; k += 1) {
const fx = BigInt(k * 4099);
const fz = BigInt(k * 6151);
assert.equal(fbm(fx, fz, SEED, 1, SHIFT, GAIN), valueNoise(fx, fz, SEED, SHIFT));
}
});
test("负例：漏掉归一化（不除以幅度和）会给出更大的值", () => {
// 故意写错：只累加 fmul(n, amp)，不做 fdiv(total, norm)
const wrong = (fx: bigint, fz: bigint): bigint => {
let total = 0n;
let amp = ONE;
for (let o = 0; o < 4; o += 1) {
total += fmul(valueNoise(fx, fz, SEED + BigInt(o), SHIFT - BigInt(o)), amp);
amp = fmul(amp, GAIN);
}
return total;
};
const fx = 12345n;
const fz = 67890n;
const good = fbm(fx, fz, SEED, 4, SHIFT, GAIN);
const bad = wrong(fx, fz);
assert.ok(bad > good, `漏归一化的值 ${bad} 未大于正确值 ${good}`);
assert.ok(bad < 2n * good, "两者差距过大，负例可能选点不当");
});
test("fBm 输出恒在 [0, ONE)", () => {
for (let k = 0; k < 300; k += 1) {
const v = fbm(BigInt(k * 2311), BigInt(k * 3517), SEED, 8, SHIFT, GAIN);
assert.ok(v >= 0n && v < ONE, `fBm 越界 ${v}`);
}
});
四条断言里最有价值的是第二条：单层 fBm 必须退化为值噪声。因为 octaves=1 时 total = fmul(n, ONE) = n、norm = ONE、fdiv(n, ONE) = n，结果应当逐位等于 valueNoise。这条断言把「归一化」与「叠加」两段代码的正确性一次性钉死——任何一处算错，单层都退不出来。第三条是负例：故意去掉归一化，值必须变大（实测确实变大且不到两倍），证明 digest 之外的判据也有区分力。
阅读导航：上一节：9.1 值噪声与 Perlin · 下一节：9.3 程序化纹理的整数化 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。