1.2 Phong 与高光
上一节只算了漫反射：光线打在粗糙表面上向各方向散开，与观察方向无关。可现实里的物体大多有一层光泽——金属、上釉的陶瓷、湿润的皮肤——它们在某个特定角度上会亮起来。这个亮点与观察方向强相关，是上一节的模型给不出的。
站内 /graphics-lighting-shadows/ 那篇把高光当作 PBR 的一部分，用微表面分布函数去描述；本卷先走经验路线，把 Bui Tuong Phong 在 1975 年给出的那个模型写成可断言的整数，看它的每一项各花多少运算，以及它和后来 Blinn 的改法究竟差在哪里。
反射向量是一个纯整数的构造
Phong 的想法是：把光 L 关于法线 N 做镜面反射得到 R，再看 R 与观察方向 V 有多接近。反射公式写出来只有一行：R = 2(N·L)N − L。
关键在它不含任何开方或除法。N 与 L 都是单位向量，N·L 是一个标量，2(N·L)N 只是把法线按这个标量放大，再减去 L。整条链是三次 fmul、三次减法，全部整数：
# snippet: 节选自 lantern/sim/cg_shade_01_2_phong.py 的 phong_spec()
t = o.dot(n, l) * 2 # 2(N·L)，一个 Q16.16 标量
r = tuple(o.m(t, n[i]) - l[i] for i in range(3)) # R = 2(N·L)N − L，纯整数
rv = max(0, o.dot(r, v)) # R·V，符号由整数比较得出
return pow_chain(o, rv, s), rv
幂次用平方链
高光的锐度由指数 s 控制：spec = max(0, R·V)^s。s 越大，亮点越小越集中。朴素写法是把 R·V 自乘 s 次，s = 32 就是 31 次乘法（连乘初值 1.0 那次共 32 次）。但定点乘法是结合律成立的精确整数运算，可以用平方链（exponentiation by squaring）：
# snippet: 节选自 lantern/sim/cg_shade_01_2_phong.py 的 pow_chain()
r, b = None, x
while e:
if e & 1:
r = b if r is None else o.m(r, b)
e >>= 1
if e:
b = o.m(b, b)
return cglib.ONE if r is None else r
s 是 2 的幂时，二进制位里只有最高位是 1，于是链上只发生自乘：x^32 是五次连续平方，5 次 fmul。实测 5 次的结果与 31 次逐次相乘逐位相同——不是「近似相同」，是同一个整数。
Blinn-Phong：把反射换成半步向量
Jim Blinn 在 1977 年给出一个性质更好的替代：不构造 R，而是构造 L 与 V 的半步向量 H = normalize(L + V)，再看 N 与 H 的夹角：
# snippet: 节选自 lantern/sim/cg_shade_01_2_phong.py 的 blinn_spec()
h = o.norm(tuple(l[i] + v[i] for i in range(3))) # 这一步要归一化：1 次 fsqrt + 3 次 fdiv
ndh = max(0, o.dot(n, h))
return pow_chain(o, ndh, s), ndh
H 的几何意义是：如果一个表面的法线恰好是 H，它就会把 L 镜面反射进 V。所以 N·H 衡量的是「真实法线离那个理想法线有多远」。这正是微表面理论后来采用的写法（第 6 章），也是 Blinn-Phong 被认为更「正确」的原因。但代价也要说清：normalize 需要一次 fsqrt 加三次 fdiv。
完整脚本：
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_01_2_phong.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_01_2_phong.py
"""1.2 Phong 与高光：反射向量的整数构造、幂次的平方链、半步向量的替代。
固定 N=(0,0,1) 与一盏方向光 L，把观察方向 V 扫过八面体细分球的 66 个顶点：
Phong R = 2(N·L)N − L spec = max(0, R·V)^s
Blinn-Phong H = normalize(L + V) spec = max(0, N·H)^s
R 是纯整数构造（无开方、无除法）；H 必须归一化（1 次 fsqrt + 3 次 fdiv），这就是两者的
全部代价差。幂次走平方链：x^32 = 5 次 fmul，与 31 次逐次相乘逐位相同。再扫一遍
shininess ∈ {4,8,16,32,64}，看叶宽如何随指数收窄，并实测一条常被说反的直觉：
低指数下 Phong 与 Blinn-Phong **都会**在 N·V<0 的方向上给出非零高光。
产物：out/cg_shade_01_2_phong.csv、.json、.svg
"""
import json
import os
import sys
import xml.etree.ElementTree as ET
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
SPEC = os.path.join(ROOT, "config", "cg_shade_01_2_phong.json")
class Ops:
"""带计数的定点运算：一次高光求值用了多少次 fmul/fdiv/fsqrt（真跑，不是估算）。"""
__slots__ = ("mul", "div", "sqrt")
def __init__(self):
self.mul = self.div = self.sqrt = 0
def m(self, a, b):
self.mul += 1
return cglib.fmul(a, b)
def d(self, a, b):
self.div += 1
return cglib.fdiv(a, b)
def r(self, x):
self.sqrt += 1
return cglib.fsqrt(x)
def dot(self, a, b):
return self.m(a[0], b[0]) + self.m(a[1], b[1]) + self.m(a[2], b[2])
def norm(self, a):
n = self.r(self.dot(a, a))
return (0, 0, 0) if n == 0 else (self.d(a[0], n), self.d(a[1], n), self.d(a[2], n))
def pow_chain(o, x, e):
"""整数平方链 x^e：e 是 2 的幂时恰好 log2(e) 次自乘。"""
r, b = None, x
while e:
if e & 1:
r = b if r is None else o.m(r, b)
e >>= 1
if e:
b = o.m(b, b)
return cglib.ONE if r is None else r
def phong_spec(o, n, l, v, s):
"""R = 2(N·L)N − L（纯整数），spec = max(0, R·V)^s。"""
t = o.dot(n, l) * 2
r = tuple(o.m(t, n[i]) - l[i] for i in range(3))
rv = max(0, o.dot(r, v))
return pow_chain(o, rv, s), rv
def blinn_spec(o, n, l, v, s):
"""H = normalize(L + V)，spec = max(0, N·H)^s。"""
h = o.norm(tuple(l[i] + v[i] for i in range(3)))
ndh = max(0, o.dot(n, h))
return pow_chain(o, ndh, s), ndh
def spec_svg(path, points):
"""横轴 N·V、纵轴高光系数的极简 SVG（纯 xml.etree，无第三方绘图库）。"""
W, H, M = 640, 400, 40
px = lambda x: M + (x + 1.0) / 2.0 * (W - 2 * M)
py = lambda y: H - M - min(1.0, max(0.0, y)) * (H - 2 * M)
svg = ET.Element("svg", {"xmlns": "http://www.w3.org/2000/svg", "width": str(W),
"height": str(H), "viewBox": f"0 0 {W} {H}"})
ET.SubElement(svg, "rect", {"x": "0", "y": "0", "width": str(W), "height": str(H),
"fill": "#0e1016"})
ET.SubElement(svg, "line", {"x1": str(M), "y1": f"{py(0):.1f}", "x2": str(W - M),
"y2": f"{py(0):.1f}", "stroke": "#2a2f3a"})
ET.SubElement(svg, "line", {"x1": f"{px(0):.1f}", "y1": str(M), "x2": f"{px(0):.1f}",
"y2": str(H - M), "stroke": "#2a2f3a"})
for key, color in (("p", "#ff6b6b"), ("b", "#5b8cff")):
seq = sorted(points, key=lambda q: q["ndotv"])
ET.SubElement(svg, "polyline", {"fill": "none", "stroke": color, "stroke-width": "2",
"points": " ".join(f"{px(q['ndotv']):.1f},"
f"{py(q[key]):.1f}" for q in seq)})
ET.ElementTree(svg).write(path, encoding="utf-8", xml_declaration=True)
def main():
spec = json.load(open(SPEC, encoding="utf-8"))
n, l = tuple(spec["normalQ16"]), tuple(spec["directionalLight"]["dirQ16"])
s = spec["shininess"]
verts, _f = cglib.build_sphere((0, 0, 0), cglib.ONE, spec["viewSphereSubdiv"])
views = [cglib.v_norm(v) for v in verts]
x = cglib.fmul(cglib.v_dot(l, l), cglib.ONE) # 任意 Q16.16 测试值
o = Ops()
chain = pow_chain(o, x, s)
chain_mul = o.mul
repeated = cglib.ONE
for _ in range(s):
repeated = cglib.fmul(repeated, x)
rows, points = [], []
for i, v in enumerate(views):
sp, _ = phong_spec(Ops(), n, l, v, s)
sb, _ = blinn_spec(Ops(), n, l, v, s)
ndotv = cglib.v_dot(n, v)
rows.append({"i": i, "v": v, "ndotv": ndotv, "sp": sp, "sb": sb})
points.append({"ndotv": ndotv / cglib.ONE, "p": sp / cglib.ONE, "b": sb / cglib.ONE})
sp_arr = np.array([r["sp"] for r in rows], dtype=np.int64)
sb_arr = np.array([r["sb"] for r in rows], dtype=np.int64)
nv_arr = np.array([r["ndotv"] for r in rows], dtype=np.int64)
diff = np.abs(sp_arr - sb_arr)
back = nv_arr < 0
sweep = []
for e in (4, 8, 16, 32, 64):
pn = bn = ob = bb = mx = 0
for v in views:
sp, _ = phong_spec(Ops(), n, l, v, e)
sb, _ = blinn_spec(Ops(), n, l, v, e)
nv = cglib.v_dot(n, v)
pn += sp > 0
bn += sb > 0
ob += nv < 0 and sp > 0
bb += nv < 0 and sb > 0
mx = max(mx, abs(sp - sb))
sweep.append({"shininess": e, "phong_nonzero": pn, "blinn_nonzero": bn,
"phong_off_specular": ob, "blinn_back_face": bb, "max_abs_diff_q16": mx})
cp, cb = Ops(), Ops()
phong_spec(cp, n, l, views[0], s)
blinn_spec(cb, n, l, views[0], s)
out = os.path.join(ROOT, "out")
os.makedirs(out, exist_ok=True)
cglib.write_csv(
os.path.join(out, "cg_shade_01_2_phong.csv"),
["index", "vx", "vy", "vz", "ndotv_q16", "ndotv", "phong_q16", "phong",
"blinn_q16", "blinn", "diff_q16"],
[[r["i"], r["v"][0], r["v"][1], r["v"][2], r["ndotv"], f"{r['ndotv'] / cglib.ONE:.6f}",
r["sp"], f"{r['sp'] / cglib.ONE:.6f}", r["sb"], f"{r['sb'] / cglib.ONE:.6f}",
abs(r["sp"] - r["sb"])] for r in rows])
spec_svg(os.path.join(out, "cg_shade_01_2_phong.svg"), points)
stats = {
"shininess": s, "normalQ16": list(n), "lightDirQ16": list(l), "views": len(rows),
"pow_chain_mul": chain_mul, "repeated_mul": s,
"pow_chain_value_q16": chain, "repeated_value_q16": repeated,
"pow_chain_equals_repeated": chain == repeated,
"phong_nonzero": int((sp_arr > 0).sum()), "blinn_nonzero": int((sb_arr > 0).sum()),
"max_abs_diff_q16": int(diff.max()), "mean_abs_diff_q16": round(float(diff.mean()), 4),
"differing_views": int((diff > 0).sum()), "back_facing_views": int(back.sum()),
"phong_off_specular": int((back & (sp_arr > 0)).sum()),
"blinn_on_back_face": int((back & (sb_arr > 0)).sum()),
"shininess_sweep": sweep,
"phong_ops": {"mul": cp.mul, "div": cp.div, "sqrt": cp.sqrt},
"blinn_ops": {"mul": cb.mul, "div": cb.div, "sqrt": cb.sqrt},
"golden": [{"v": list(r["v"]), "sp": r["sp"], "sb": r["sb"], "ndotv": r["ndotv"]}
for r in rows],
}
with open(os.path.join(out, "cg_shade_01_2_phong.json"), "w", encoding="utf-8") as fh:
json.dump(stats, fh, ensure_ascii=False, indent=2)
print(f"观察方向数 : {stats['views']} 背面方向数 {stats['back_facing_views']}")
print(f"平方链 x^{s} : {chain_mul} 次 fmul vs 逐次相乘 {s} 次，"
f"相同 ? {stats['pow_chain_equals_repeated']}")
print(f"高光非零方向数 : Phong {stats['phong_nonzero']} / "
f"Blinn-Phong {stats['blinn_nonzero']}")
print(f"两者最大差 : {stats['max_abs_diff_q16']} Q16.16 平均差 "
f"{stats['mean_abs_diff_q16']} 不同方向数 {stats['differing_views']}")
print(f"Phong 离镜面 : {stats['phong_off_specular']} "
f"Blinn 背面高光 : {stats['blinn_on_back_face']}")
for r in sweep:
print(f" s={r['shininess']:3d} phong={r['phong_nonzero']:3d} "
f"blinn={r['blinn_nonzero']:3d} off_spec={r['phong_off_specular']:2d} "
f"blinn_back={r['blinn_back_face']:2d} maxdiff={r['max_abs_diff_q16']}")
print(f"一次求值运算量 : Phong mul={cp.mul} div={cp.div} sqrt={cp.sqrt} | "
f"Blinn mul={cb.mul} div={cb.div} sqrt={cb.sqrt}")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一次：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_shade_01_2_phong.py
观察方向数 : 66 背面方向数 25
平方链 x^32 : 5 次 fmul vs 逐次相乘 32 次，相同 ? True
高光非零方向数 : Phong 12 / Blinn-Phong 18
两者最大差 : 31240 Q16.16 平均差 1190.8333 不同方向数 18
Phong 离镜面 : 0 Blinn 背面高光 : 0
s= 4 phong= 27 blinn= 53 off_spec= 4 blinn_back=12 maxdiff=33480
s= 8 phong= 27 blinn= 45 off_spec= 4 blinn_back= 4 maxdiff=31665
s= 16 phong= 15 blinn= 31 off_spec= 0 blinn_back= 2 maxdiff=27804
s= 32 phong= 12 blinn= 18 off_spec= 0 blinn_back= 0 maxdiff=31240
s= 64 phong= 3 blinn= 9 off_spec= 0 blinn_back= 0 maxdiff=24370
一次求值运算量 : Phong mul=14 div=0 sqrt=0 | Blinn mul=11 div=3 sqrt=1
三组数字
平方链 5 次 fmul，逐次相乘 32 次，结果相同。 省下 27 次乘法，而且省得随指数增长：s = 64 只需 6 次。
叶宽：Phong 12 个方向非零，Blinn-Phong 18 个。 同一 shininess 下 Blinn-Phong 的叶更宽。原因在 H 的定义里：N·H > 0 等价于 N·V > −N·L，这个门槛比 R·V > 0 宽松，于是高光多盖住一圈方向。
一次求值的运算量：Phong mul=14 div=0 sqrt=0，Blinn mul=11 div=3 sqrt=1。 这组数字推翻一句流行话——「Blinn-Phong 更省」。它确实少 3 次乘法，但换来 1 次平方根和 3 次除法。在定点内核里，fmul 是一次 64 位整数乘加一次右移，fdiv 是 64 位整数除法，fsqrt 是牛顿迭代的 isqrt——后两者的单次代价都远高于乘法。所以在本工程的口径下，逐像素求值时 Blinn-Phong 反而更贵。（若把 H 放到顶点阶段算、再插值到像素，结论会反过来；这正是第 13 章「着色变体」要讨论的取舍。）
一条被说反的直觉
常有人说「Blinn-Phong 不会照亮背面」。扫描结果说不是：
s= 4 phong= 27 blinn= 53 off_spec= 4 blinn_back=12 maxdiff=33480
s= 8 phong= 27 blinn= 45 off_spec= 4 blinn_back= 4 maxdiff=31665
s= 16 phong= 15 blinn= 31 off_spec= 0 blinn_back= 2 maxdiff=27804
off_spec 是 N·V < 0（观察者在背面）却 spec > 0 的方向数。s = 4 时 Phong 有 4 个、Blinn-Phong 有 12 个——两个都会，而且低指数下 Blinn 更多。指数升到 32 以上，两者的背面高光才都消失（off_spec = 0、blinn_back = 0）。
真正的区别在别处：Blinn-Phong 的叶更宽，掠射角附近更早出现高光；Phong 的叶更窄更尖，在低模上容易出现「跳变」。至于「谁更正确」，要到第 6 章的微表面理论才能回答——那里会把 N·H 换成分布函数 D，把 max(0, N·L) 换成几何项与菲涅尔项。
跨语言核对
P 层用 BigInt 重算 66 个方向的两种高光，与 S 侧 golden 逐项比对，并复现叶宽与最大差：
// snippet: 节选自 lantern/core/test/cg_shade_01_2_phong.test.ts
const x = fmul(vDot(l, l), ONE);
let repeated = ONE;
for (let i = 0; i < s; i += 1) repeated = fmul(repeated, x);
assert.equal(powChain(x, s), repeated); // 平方链 == 逐次相乘，逐位相同
// 负例：反射向量写成 R = 2(N·L)N + L，高光必须与 golden 分叉
阅读导航：上一节：1.1 Lambert 与漫反射 · 下一节：1.3 着色频率与插值 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。