本文合并归档了 2020 年保存的《用 C 语言画光》六篇笔记，按原顺序保留正文、代码与插图，并整理标题和 Markdown 格式。原系列及示例代码来源：miloyip/light2d。文中的第一人称叙述沿用原笔记，并非本站作者的原创经历。
用 C 语言画光（一）：基础
1. 光
假设我们在一个二维的世界，这里有些会发光的二维形状，并暂时只考虑单色光。我们想知道的是，在这个空间中，每一点从 360 度各方向共有多少光经过。换成数学方式表示，我们想对每个二维坐标 $(x,y)$ 求积分：
$$ F(x,y) = \int_{0}^{2\pi}L(x,y,\theta)d\theta $$当中 $L(x,y,\theta)$ 代表在二维坐标 $(x,y)$ 在 $\theta$ 方向有多少光经过。我们想采样这个 $F(x,y)$ 存储在图像中，并利用 svpng() 输出结果，那么我们的 main() 函数简单为：
#include "svpng.inc"
#include <math.h>
#define W 512
#define H 512
unsigned char img[W * H * 3];
// ...
int main() {
unsigned char* p = img;
for (int y = 0; y < H; y++)
for (int x = 0; x < W; x++, p += 3)
p[0] = p[1] = p[2] = (int)(fminf(sample(
(float)x / W, (float)y / H) * 255.0f, 255.0f));
svpng(fopen("basic.png", "wb"), W, H, img, 0);
}
无论图像长宽是多少，这个二维空间的坐标范围都是 $(x,y) \in [0, 1] \times [0, 1]$ 。我们把结果映射至 $\{0, 1,\ldots,255\}$ 。
2. 蒙地卡罗积分
由于无法以解析式求解这个积分，我们使用蒙地卡罗积分法（Monte Carlo integration）。
在这个问题中，我们随机采样 $N$ 个方向 $\theta_1,\theta_2,\ldots,\theta_N$ ，然后计算 $L(x,y,\theta_i)$ 的平均值：
代码实现也很浅白（设每像素向 64 个随机方向均匀采样／uniform sampling）：
#define TWO_PI 6.28318530718f
#define N 64
float sample(float x, float y) {
float sum = 0.0f;
for (int i = 0; i < N; i++) {
float a = TWO_PI * rand() / RAND_MAX;
sum += trace(x, y, cosf(a), sinf(a));
}
return sum / N;
}
当中 trace(ox, oy, dx, dy) 函数代表从 $o$ 位置从单位矢量 $\hat{d}$ 方向接收到的光。
（更新：本文不考虑实际单位，所以实现时把系数 $2\pi$ 去掉了。感谢 @Bimos 提醒）
3. 光线步进
通常，我们可以用光线追踪（ray tracing）方法，求出光线 $r(t)=o+\hat{d}t$ 与场景的最近点。
然而，我们需要为每种几何形状编写与光线的相交函数（通常比较复杂），之后做一些效果可能还要提供相交点的法线（normal vector）。
为简单起见，本文采用光线步进（ray marching）方法（又称为球体追踪／sphere tracing [1]），场景只需要以带符号距离场（signed distance field, SDF） 表示
- 当 $\phi(x) \gt 0$ ，表示坐标 $x$ 位于场景形状之外，且 $x$ 与最近形状边界的距离为 $\phi(x)$ ；
- 当 $\phi(x) \lt 0$ ，表示坐标 $x$ 位于场景形状之内，且 $x$ 与最近形状边界的距离为 $-\phi(x)$ ；
- 当 $\phi(x) = 0$ ，说明 $x$ 刚好在形状边界上。
例如，圆心为 $c$ 、半径为 $r$ 的圆形 SDF 定义为（更精确的说法是圆盘／disk）：
$$ \phi_{\text{circle}}(x)=\|x-c\|-r $$用代码表示：
float circleSDF(float x, float y, float cx, float cy, float r) {
float ux = x - cx, uy = y - cy;
return sqrtf(ux * ux + uy * uy) - r;
}
而所谓光线步进，就是我们从光线起始点 $t=0$ ，逐步增加 $t$ ，当 $\phi(r(t)) \le 0$ 代表我们到达到场景中某个形状的表面或内部。那么我们每次可以步进多远？由于 $\phi(x) \gt 0$ 时，代表 $x$ 距最近形状的距离，所以我们至少可以步进该距离而不会碰到任何形状！
设场景只有一个发光的圆形，圆心为 $(0.5, 0.5)$ ，半径为 $0.1$ ，它每点都向各方向发射 $2$ 个单位的光。那么光线步进可实现为：
#define MAX_STEP 10
#define MAX_DISTANCE 2.0f
#define EPSILON 1e-6f
float trace(float ox, float oy, float dx, float dy) {
float t = 0.0f;
for (int i = 0; i < MAX_STEP && t < MAX_DISTANCE; i++) {
float sd = circleSDF(ox + dx * t, oy + dy * t, 0.5f, 0.5f, 0.1f);
if (sd < EPSILON)
return 2.0f;
t += sd;
}
return 0.0f;
}
如果光线超过指定距离，或是步数太多，都终止步进。注意我们只能尽量接近表面，所以用 表示足够近的阈值。整个程序完成，运算结果如下：
可以看到图像中有很多噪点，这是由于蒙地卡罗积分法具随机性，计算出来的估值与精确数值会有误差。增加采样数目 $N$ ，就能令结果更精确（准确地说是减少方差／variance）：
然而，随着 $N$ 上升，运算时间也线性上升。有没有方法可以改善？
4. 分层、抖动采样
形成噪点的原因在于，每个像素所得到的随机方向都很不一样。那么，如果我们不用随机方向，而是平分 360 度向 $N$ 个方向采样，效果会如何？这种方式称为分层采样（stratified sampling）：
float sample(float x, float y) {
float sum = 0.0f;
for (int i = 0; i < N; i++) {
// float a = TWO_PI * rand() / RAND_MAX; // 均匀采样
float a = TWO_PI * i / N; // 分层采样
// ...
改变这一行代码的结果是：
很好，没有噪点，但也太过规律了！
我们可以结合上面两种采样方式，就是先分为 $N$ 个角度区域，再从区域中均匀采样，这称为抖动采样（jittered sampling）。
// float a = TWO_PI * rand() / RAND_MAX; // 均匀采样
// float a = TWO_PI * i / N; // 分层采样
float a = TWO_PI * (i + (float)rand() / RAND_MAX) / N; // 抖动采样
改变这一行代码的结果是：
同样使用每像素 $N=64$ 次采样，仅仅改变采样方式（一行代码），效果就好很多：
5. 结语
这个完整程序位于 basic.c。如果不含加载头文件及常数定义，仅有 30 行代码。你可以改一下圆形位置、大小、光度，测试时可用较小的 W、H 和 N 加速运行过程。
本文简单介绍了蒙地卡罗积分、光线步进、分层采样等概念。下一篇《构造实体几何》会讲如何在场景中加入更多形状，将会显示出阴影。而之后也会尝试实现一些光学法则，生成更有趣的图形。
参考
[1] Hart, John C. “Sphere tracing: A geometric method for the antialiased ray tracing of implicit surfaces.” The Visual Computer12.10 (1996): 527-545.
用 C 语言画光（二）：构造实体几何
在《用 C 语言画光（一）：基础》中，我们画了一个发光的圆形，那么，该如何画更多形状呢？
1. 场景函数
重温一下，之前的光线步进实现是这样的：
float trace(float ox, float oy, float dx, float dy) {
float t = 0.0f;
for (int i = 0; i < MAX_STEP && t < MAX_DISTANCE; i++) {
float sd = circleSDF(ox + dx * t, oy + dy * t, /*...*/); // SDF
if (sd < EPSILON)
return 2.0f; // 自发光
t += sd;
}
return 0.0f;
}
在步进的过程中，我们要取得当前采样点与场景形状的带符号距离（signed distance），当少于阈值，就要返回该形状的自发光强度（emissive）。我们首先进行代码重构（code refactoring），把场景采样功能提取成一个函数（extract function），它返回 SDF 和自发光。
typedef struct { float sd, emissive; } Result;
Result scene(float x, float y) {
Result r = { circleSDF(x, y, 0.5f, 0.5f, 0.1f), 2.0f };
return r;
}
float trace(float ox, float oy, float dx, float dy) {
float t = 0.0f;
for (int i = 0; i < MAX_STEP && t < MAX_DISTANCE; i++) {
Result r = scene(ox + dx * t, oy + dy * t); // <-
if (r.sd < EPSILON) // <-
return r.emissive; // <-
t += r.sd; // <-
}
return 0.0f;
}
那么，我们之后改动场景时就不用更改 trace() 函数，只需修改 scene() 。
2. 加入多个形状
我们在上一篇谈到，$\phi(x)$ 表示坐标 $x$ 与场景形状的最近距离。如果场景中有两个形状 $A$ 和 $B$，它们的 SDF 分别为为 $\phi_A(x)$ 和 $\phi_B(x)$ ，那么坐标 $x$ 与这两个形状的最近距离，就是两个 SDF 的最小值，即 $\min(\phi_A(x), \phi_B(x))$。这个操作可以定义为，形状 $A$ 和 $B$ 的并集（union）$A \cup B$，其 SDF 为：
$$ \phi_{A\cup B}(x) = \min(\phi_A(x), \phi_B(x)) \\ \phi_{A\cap B}(x) = \max(\phi_A(x), \phi_B(x)) \\ \phi_{A \backslash B}(x) = \max(\phi_A(x), -\phi_B(x)) \\ $$实现也很简单：
Result unionOp(Result a, Result b) {
return a.sd < b.sd ? a : b;
}
我们加入 3 个圆形试试，每个圆形的自发光值都不同：
Result scene(float x, float y) {
Result r1 = { circleSDF(x, y, 0.3f, 0.3f, 0.10f), 2.0f };
Result r2 = { circleSDF(x, y, 0.3f, 0.7f, 0.05f), 0.8f };
Result r3 = { circleSDF(x, y, 0.7f, 0.5f, 0.10f), 0.0f };
return unionOp(unionOp(r1, r2), r3);
}
运行结果：
我们发现，此图形增添了一个视觉上的特徵──阴影！这是因为，进行光线步进时，追踪至非常接近形状表面时便会停下来，不会继续穿透形状，所以形状之间是有遮挡关系的，也就形成了阴影。而且因为光源是线条（此例子中是圆孤边界）而不是一点，所以会产生半影（penumbra）的效果。
3. 构造实体几何
采用 SDF 来表示形状，还可以简单地实现构造实体几何（constructive solid geometry），即是以形状点集的布尔操作来表示模型。3 个基本运算为并集（union）、交集（intersection）和相对补集（relative complement，或称为 set difference）：
注意，相对补集运算为非可交换的（non-commutative）。另外，涉及凹的形状时，得出的 SDF 并不代表最近距离，但它的估值是保守的，我们仍然可以用光线步进方法来找出边界。
实现时，要注意除了计算结果的 SDF，还要考虑采用哪个形状的特性（暂时只有自发光）。对于交集，我们要采用相反的形状；对于相对补集则采用左运算元的形状：
Result unionOp(Result a, Result b) {
return a.sd < b.sd ? a : b;
}
Result intersectOp(Result a, Result b) {
Result r = a.sd > b.sd ? b : a;
r.sd = a.sd > b.sd ? a.sd : b.sd;
return r;
}
Result subtractOp(Result a, Result b) {
Result r = a;
r.sd = (a.sd > -b.sd) ? a.sd : -b.sd;
return r;
}
做一个测试场景，左方是较亮的圆形，右方是较暗的圆形：
Result scene(float x, float y) {
Result a = { circleSDF(x, y, 0.4f, 0.5f, 0.20f), 1.0f };
Result b = { circleSDF(x, y, 0.6f, 0.5f, 0.20f), 0.8f };
return unionOp(a, b);
// return intersectOp(a, b);
// return subtractOp(a, b);
// return subtractOp(b, a);
}
把 4 个语句逐一取消注释，运行结果：
4. 结语
通过 CSG，可以把各种形状组合成不同的形状，这在运算在工业用的建模工具非常常见。你可以尝试做不同的形状，例如齿轮、花朵等。下一篇会讲述更多形状的 SDF 定义和相关运算。
本文的完整代码位于 csg.c。
用 C 语言画光（三）：形状
再前两篇中，我们一直只用到圆形作为基础形状，本文再介绍一些常用形状的 SDF。但这些形状也只是作为例子，重要的是掌握如何推导出这些 SDF。
1. 圆形
我们先重温一下圆形的 SDF：
$$ \phi_{\text{circle}}(x)=\|x-c\|-r $$当中 $\|x-c\|$ 就是坐标 $x$ 与一个点 $c$（圆心）的欧氏距离（Euclidean distance）。把这个距离减去 $r$，实际上等同于把一点向外扩展 $r$ 的范围。其实，对于任何形状的 SDF，减去一个常数 $r$ ，都等同于把该形状它向外扩展 $r$ 的范围。本文稍后也会看到更多应用例子。
在数学上，我们还可考虑这是任何形状与（圆心位于原点的）圆盘的闵可夫斯基和（Minkowski sum）。
2. 平面
除了圆形，最简单的形状应该是平面（plane）。在二维中这个形状应该称作直线，或者更精确地应称作二维的半空间（half-space），就是把二维空间用直线切割成两个无限大面积的部分。我们这里暂且称这个形状为平面吧。
任意维度空间的（超）平面都可以用平面上的一点 $p$ 与法线（normal vector）$n$ 来定义，由于平面与法线垂直，平面上的任意一个矢量 $x-p$ 与法线的点积为零：
$$ (x-p) \cdot n = 0 $$如果法矢量为单位矢量 $\hat{n}$，空间中任意一点与平面的带符号距离为：
$$ \phi_{\text{plane}}(x) = (x-p) \cdot \hat{n} = 0 $$我们还可把点积展开成这种形式：
$$ \phi_{\text{plane}}(x) = Ax + By +c $$但为了方便使用，我们的实现用点和法线作为参数：
float planeSDF(float x, float y, float px, float py, float nx, float ny) {
return (x - px) * nx + (y - py) * ny;
}
Result scene(float x, float y) {
Result b = { planeSDF(x, y, 0.0f, 0.5f, 0.0f, 1.0f), 0.8f };
return b;
}
这个例子中，上半部分是发射 0.8 自发光的形状：
我们可以用 CSG，例如把圆形与上面的平面相交，得出半圆：
Result scene(float x, float y) {
Result a = { circleSDF(x, y, 0.5f, 0.5f, 0.2f), 1.0f };
Result b = { planeSDF(x, y, 0.0f, 0.5f, 0.0f, 1.0f), 0.8f };
return intersectOp(a, b);
}
3. 胶囊形
在《用 C 语言画直线》已使用到胶囊形（capsule）的 SDF，在这里简单推导一下。
我们把胶囊形以两点 $\mathbf{a}, \mathbf{b}$ 和半径 $r$ 来表示：
有些场合可能每使用一点 $\mathbf{a}$ 及矢量 $\mathbf{u} = \mathbf{b} - \mathbf{a}$ 去表示，实际上我们在计算其 SDF 时也需要做这个减法。
计算胶囊形的 SDF，等价于计算一点与线段 ab 的距离，再扩大半径 $r$（即把结果减去 $r$ ）。而计算一点与线段的距离的方法，是把点 $\mathbf{x}$ 投影在线段的延长线上，然后把投影点 $\mathbf{x}^\prime$ 控制在线段的范围内：
上式把最后的 $\frac{1}{\|\mathbf{u}\|}$ 因子放进括号里，减少开方运算。然后，所需的距离为：
$$ d = \|\mathbf{p} - \mathbf{p}^\prime\| $$这里我们把线段和胶囊形的 SDF 分开为两个函数实现：
float segmentSDF(float x, float y, float ax, float ay, float bx, float by) {
float vx = x - ax, vy = y - ay, ux = bx - ax, uy = by - ay;
float t = fmaxf(fminf((vx * ux + vy * uy) / (ux * ux + uy * uy), 1.0f), 0.0f);
float dx = vx - ux * t, dy = vy - uy * t;
return sqrtf(dx * dx + dy * dy);
}
float capsuleSDF(float x, float y, float ax, float ay, float bx, float by, float r) {
return segmentSDF(x, y, ax, ay, bx, by) - r;
}
Result scene(float x, float y) {
Result c = { capsuleSDF(x, y, 0.4f, 0.4f, 0.6f, 0.6f, 0.1f), 1.0f };
return c;
}
结果：
4. 矩形
这里实现的矩形是二维定向矩形包围盒（oriented bounding box, OBB），可以由中心点 $\mathbf{c}$、旋转角 $\theta$ 和半长（half length) $s$ 来表示：
首先我们要把点 $\mathbf{c}$ 变换至 OBB 的坐标系。我们可反过来思考，如果从 OBB 的坐标系的一点，转换至世界坐标系：
$$ \mathbf{x} = \left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta\end{array}\right] \mathbf{x}^{\prime}+\mathbf{c} $$那么，其逆变换就是（旋转矩阵是正交矩阵，有 $Q^{-1} = Q^T$ ）：
$$ \mathbf{x}^{\prime}=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\end{array}\right](\mathbf{x}-\mathbf{c}) $$因为矩形的对称性，我们只需考虑第一象限，可用绝对值把 $\mathbf{x}^{\prime}$ 移到第一象限。然后，要对内和外分计算距离，再取得最少值。
float boxSDF(float x, float y, float cx, float cy, float theta, float sx, float sy) {
float costheta = cosf(theta), sintheta = sinf(theta);
float dx = fabs((x - cx) * costheta + (y - cy) * sintheta) - sx;
float dy = fabs((y - cy) * costheta - (x - cx) * sintheta) - sy;
float ax = fmaxf(dx, 0.0f), ay = fmaxf(dy, 0.0f);
return fminf(fmaxf(dx, dy), 0.0f) + sqrtf(ax * ax + ay * ay);
}
Result scene(float x, float y) {
Result d = { boxSDF(x, y, 0.5f, 0.5f, TWO_PI / 16.0f, 0.3f, 0.1f), 1.0f };
return d;
}
结果：
减去一个半径（如 0.1）就获得圆角长方形：
Result e = { boxSDF(x, y, 0.5f, 0.5f, TWO_PI / 16.0f, 0.3f, 0.1f) - 0.1f, 1.0f };
return e;
5. 三角形
三角形的 SDF 需计算点至三个线段的最短距离，然后再判断点在内还是在外。这里我们用缠绕顺序（winding order）来定义三角形的内外。详细的判断方法请参考 [1]。
float triangleSDF(float x, float y, float ax, float ay, float bx, float by, float cx, float cy) {
float d = fminf(fminf(
segmentSDF(x, y, ax, ay, bx, by),
segmentSDF(x, y, bx, by, cx, cy)),
segmentSDF(x, y, cx, cy, ax, ay));
return (bx - ax) * (y - ay) > (by - ay) * (x - ax) &&
(cx - bx) * (y - by) > (cy - by) * (x - bx) &&
(ax - cx) * (y - cy) > (ay - cy) * (x - cx) ? -d : d;
}
Result scene(float x, float y) {
Result f = { triangleSDF(x, y, 0.5f, 0.2f, 0.8f, 0.8f, 0.3f, 0.6f), 1.0f };
return f;
}
结果：
减去一个半径就获得圆角三角形：
6. 结语
本文介绍了 5 种二维形状的 SDF，当中用到不同的技巧，例如运用形状的对称性质。长方形实际上也可用三角形的做法实现，然而不考虑对称性，所需的运算会更多。
读者可利用这些形状及 CSG 来绘画各种形状的光，也可以尝试推导其他形状的 SDF。
我们之后会再讨论用 SDF 的其他技巧，但接下来的两篇将先回到光学相关内容──反射与折射。
本文的源程序位于 shapes.c 及 m.c。
参考
[1] Dan Sunday, Inclusion of a Point in a Polygon.
用 C 语言画光（四）：反射
上一篇我们讨论了一些形状，本篇回到基于光学的模拟，先谈较简单的──光的反射（reflection）。我们假设形状的边界是平滑的，会产生镜面反射（specular reflection）。实现反射需要一点几何上的推导，以及考虑如何用 SDF 获得边界法线，最后把追踪过程改为递归方式进行。
1. 反射矢量
给定一个入射方向 $\mathbf{i}$ 及边界法线 $\hat{\mathbf{n}}$ ，我们可以根据下图计算出反射方向 $\mathbf{r}$ ：
首先，把 $\mathbf{i}$ 用点积投影至法线方向，得到法线方向的长度（注意法线需为单位矢量），然后再把这长度往法线方向缩放，得到图中蓝色虚线矢量 $(\mathbf{i} \cdot \hat{\mathbf{n}})\hat{\mathbf{n}}$ 。之后，把 $\mathbf{i}$ 减去两个蓝色虚线矢量就能得到反射矢量 $\mathbf{r}$ ：
$$ \mathbf{r} = \mathbf{i}-2(\mathbf{i} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}} $$用 C 语言实现：
void reflect(float ix, float iy, float nx, float ny, float* rx, float* ry) {
float idotn2 = (ix * nx + iy * ny) * 2.0f;
*rx = ix - idotn2 * nx;
*ry = iy - idotn2 * ny;
}
但问题是，我们怎样从场景中取得边界法线？
2. 梯度与法线
在此系列文章中，我们采用了 SDF 来构建场景。在形状的边界，距离场变化最大的方向便是法线方向。根据矢量微积分（vector calculus），一个纯量场（scalar field）的最大变化方向就是其梯度（gradient），所以我们需要求形状边界位置的 SDF 梯度：
$$ \nabla \phi=\left(\frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y}\right) $$这梯度其实就是 SDF 分别在 $x$ 和 $y$ 方向的变化率。SDF 几乎处处可导，而且，SDF 为程函方程的（Eikonal equation）特例，具特性：
$$ \|\nabla \phi\| = 1 $$换句话说，SDF 的梯度是单位矢量。
假设我们不知道 SDF 梯度的分析解，可利用数值方法，如中心差分法（central difference），求出梯度的近似值：
$$ \begin{aligned} \nabla \phi &= \left(\frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y}\right) \\ & \approx\left(\frac{\phi(x+h, y)-\phi(x-h, y)}{2 h}, \frac{\phi(x, y+h)-\phi(x, y-h)}{2 h}\right) \end{aligned} $$用 C 语言实现：
void gradient(float x, float y, float* nx, float* ny) {
*nx = (scene(x + EPSILON, y).sd - scene(x - EPSILON, y).sd) * (0.5f / EPSILON);
*ny = (scene(x, y + EPSILON).sd - scene(x, y - EPSILON).sd) * (0.5f / EPSILON);
}
我们用一个测试场景去看看它的梯度，当中有一个圆形、两个矩形：
Result scene(float x, float y) {
Result a = { circleSDF(x, y, 0.4f, 0.2f, 0.1f), 2.0f };
Result b = { boxSDF(x, y, 0.5f, 0.8f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f };
Result c = { boxSDF(x, y, 0.8f, 0.5f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f };
return unionOp(unionOp(a, b), c);
}
我们改一下 main() 函数，用颜色把梯度可视化，红色通道显示 $\frac{\partial{\phi}}{\partial{x}}$ ，绿色通道显示 $\frac{\partial{\phi}}{\partial{y}}$ ：
for (int y = 0; y < H; y++)
for (int x = 0; x < W; x++, p += 3) {
float nx, ny;
gradient((float)x / W, (float)y / H, &nx, &ny);
p[0] = (int)((fmaxf(fminf(nx, 1.0f), -1.0f) * 0.5f + 0.5f) * 255.0f);
p[1] = (int)((fmaxf(fminf(ny, 1.0f), -1.0f) * 0.5f + 0.5f) * 255.0f);
p[2] = 0;
}
上图中可以看到，大部分位置的梯度是连续的，至少我们现在需要在形状边界上的位置是连续的。
（注：由于有数值误差，计算出来的梯度分量可能略超出 $[-1, 1]$ 的范围，所以在输出时限制其范围。但这个轻微的误差不影响我们计算反射。）
3. 递归追踪
我们考虑形状边界的反射率（reflectivity）为 $[-1, 1]$ 范围，为此增加一个字段：
typedef struct { float sd, emissive, reflectivity; } Result;
Result scene(float x, float y) {
Result a = { circleSDF(x, y, 0.4f, 0.2f, 0.1f), 2.0f, 0.0f };
Result b = { boxSDF(x, y, 0.5f, 0.8f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f, 0.9f };
Result c = { boxSDF(x, y, 0.8f, 0.5f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f, 0.9f };
return unionOp(unionOp(a, b), c);
}
如果反射率超过 1，总能量就变多，不符合能量守恒。少于 1 代表形状吸收了能量。
之后，修改 trace() ，当我们在向一个方向追踪光线时，若碰到会反射的形状时，就计算该边界相交点的法线，然后递归追踪反射方向，把收到的所有光加总。为免无限反射做成无限递归，需要限制反射次数：
#define BIAS 1e-4f
#define MAX_DEPTH 3
float trace(float ox, float oy, float dx, float dy, int depth) {
float t = 0.0f;
for (int i = 0; i < MAX_STEP && t < MAX_DISTANCE; i++) {
float x = ox + dx * t, y = oy + dy * t;
Result r = scene(x, y);
if (r.sd < EPSILON) {
float sum = r.emissive;
if (depth < MAX_DEPTH && r.reflectivity > 0.0f) {
float nx, ny, rx, ry;
gradient(x, y, &nx, &ny);
reflect(dx, dy, nx, ny, &rx, &ry);
sum += r.reflectivity * trace(x + nx * BIAS, y + ny * BIAS, rx, ry, depth + 1);
}
return sum;
}
t += r.sd;
}
return 0.0f;
}
float sample(float x, float y) {
float sum = 0.0f;
for (int i = 0; i < N; i++) {
float a = TWO_PI * (i + (float)rand() / RAND_MAX) / N;
sum += trace(x, y, cosf(a), sinf(a), 0); // 开始时深度为 0
}
return sum / N;
}
另外，由于我们要向反射方向追踪，如果用原来的相交点，追踪时就为立即遇到 r.sd < EPSILON 而停止，所以我们稍微把相交点往法线方向偏移 BIAS 的距离，只要 BIAS > EPSILON ，就可以避免这个问题。但太大的话也会造成误差。结果：
我们可以看到两个矩形之间有第二次的反射。
利用平面和圆形的相对补集，就能做出凹面镜（concave mirror），弧形会产生题图中的焦散（caustics）现象：
Result scene(float x, float y) {
Result a = { circleSDF(x, y, 0.4f, 0.2f, 0.1f), 2.0f, 0.0f };
Result d = { planeSDF(x, y, 0.0f, 0.5f, 0.0f, -1.0f), 0.0f, 0.9f };
Result e = { circleSDF(x, y, 0.5f, 0.5f, 0.4f), 0.0f, 0.9f };
return unionOp(a, subtractOp(d, e));
}
4. 结语
也许大家在中学物理中都学过镜面反射的知识，但用简单的程序实现出来总会令人觉得更有趣。
本文的重点是用数值方法求 SDF 在边界上的法线，以及通过递归来实现多次反射。
下一篇将讨论折射，那么，我们就可以画出穿过透镜的光了。
本文的源代码位于 reflection.c。
用 C 语言画光（五）：折射
与上篇反射（reflection）成对出现的光学现象就是折射（refraction）。当光从一个介质（medium）进入另一个介质，便会出现这两种现象。
基于能量守恒，当没有能量散失时，反射光的能量与折射光的能量之和等于入射光的能量。在本篇中我们先假设反射和折射的比率是恒定的，例如一个形状的表面总是会反射 20% 的入射光，并折射 80% 的入射光。
但我们知道，折射的方向会改变（所以才称为「折」射），那么怎样计算折射方向呢？
1. 斯涅尔定律
读者可能在中学阶段已学会斯涅尔定律（Snell’s law），它描述了入射角、折射角与介质折射率（refractive index）的关系：
当 $\eta_1 > \eta_2$ 时，例如光线从玻璃（ $\eta_1 \approx 1.5$ ）进入空气（ $\eta_2 \approx 1$ ），会有机会出现内全反射（total internal reflection），从斯涅尔定律分析：
$$ \sin \theta_2 = \frac{\eta_1}{\eta_2} \sin \theta_1 $$如果 $\eta_1 > \eta_2 \Leftrightarrow \frac{\eta_1}{\eta_2} \gt 1$ ，某些入射角 $\theta_1$ 范围会令右侧大于 1，左侧便会超出正弦的值域 $[-1, 1]$ 。在这种情况下，不会发生折射，所有入射光会完全反射。我们把上式左侧设为 1，则可求出临界角（critical angle）：
$$ \begin{aligned} &1=\frac{\eta_{1}}{\eta_{2}} \sin \theta_{c}\\ &\Leftrightarrow \theta_{c}=\sin ^{-1} \frac{\eta_{2}}{\eta_{1}} \end{aligned} $$2. 折射矢量
我们以矢量表示光线方向，需要求出折射方向。这会比上一篇反射复杂一点点。先看图：
与反射不一样，折射最好使用单位矢量表示入射方向 $\hat{\mathbf{i}}$ 和折射方向 $\hat{\mathbf{t}}$ 。这里还定义了界面的切线方向 $\hat{\mathbf{m}}$ 。从上图可以看到，折射矢量可从法线方向与切线方向的矢量之和得出：
$$ \hat{\mathbf{t}} = \sin \theta_2 \hat{\mathbf{m}} - \cos\theta_2 \hat{\mathbf{n}} $$接下来，我们要把上式各部分以 $\hat{\mathbf{i}},\hat{\mathbf{n}}, \eta_1, \eta_2$ 表示。切线 $\hat{\mathbf{m}}$ 可以从入射方向 $\hat{\mathbf{i}}$ 分解得出：
$$ \begin{aligned} \hat{\mathbf{i}} &=(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}+\sin \theta_{1} \hat{\mathbf{m}} \\ \Leftrightarrow \hat{\mathbf{m}} &=\frac{\hat{\mathbf{i}}-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}}{\sin \theta_{1}} \end{aligned} $$然后 $\sin \theta_2$ 可通过斯涅尔定律以 $\sin\theta_1, \eta_1, \eta_2$ 表示，而 $\sin\theta_1$ 又可以转换成点积：
$$ \sin\theta_2 = \frac{\eta_1}{\eta_2}\sin\theta_1=\frac{\eta_1}{\eta_2}\sqrt{1-{\cos^2\theta_1}} = \frac{\eta_1}{\eta_2}\sqrt{1-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}})^2} $$有正弦也就可以求余弦：
$$ \cos\theta_2 = \sqrt{1 - \sin^2\theta_2} = \sqrt{1 - {\left(\frac{\eta_1}{\eta_2}\right)}^2 \left(1-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}})^2\right)} $$把(2)、(3)、(4) 代入 (1)，简化：
$$
\begin{aligned} \hat{\mathbf{t}} &=\sin \theta_{2} \hat{\mathbf{m}}-\cos \theta_{2} \hat{\mathbf{n}} \ &=\frac{\sin \theta_{2}}{\sin \theta_{1}}(\hat{\mathbf{i}}-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}})-\sqrt{1-\left(\frac{\eta_{1}}{\eta_{2}}\right)^{2}\left(1-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}})^{2}\right)} \hat{\mathbf{n}} \ &=\frac{\eta_{1}}{\eta_{2}}(\hat{\mathbf{i}}-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}})-\sqrt{1-\left(\frac{\eta_{1}}{\eta_{2}}\right)^{2}\left(1-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}})^{2}\right)} \hat{\mathbf{n}} \ &=\frac{\eta_{1}}{\eta_{2}} \hat{\mathbf{i}}-\left(\frac{\eta_{1}}{\eta_{2}}(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}})+\sqrt{1-\left(\frac{\eta_{1}}{\eta_{2}}\right)^{2}\left(1-(\hat{\mathbf{i}} \cdot \hat{\mathbf{n}})^{2}\right)}\right) \hat{\mathbf{n}} \end{aligned}
$$
噢，计算折射矢量竟然一个三角函数都不需要！由于上式只需折射率之比值，我们的函数使用 $\eta=\eta_1 / \eta_2$ 作为参数，实现如下：
int refract(float ix, float iy, float nx, float ny, float eta, float* rx, float* ry) {
float idotn = ix * nx + iy * ny;
float k = 1.0f - eta * eta * (1.0f - idotn * idotn);
if (k < 0.0f)
return 0; // 全内反射
float a = eta * idotn + sqrtf(k);
*rx = eta * ix - a * nx;
*ry = eta * iy - a * ny;
return 1;
}
当根号内的值为负数，也是代表全内反射，我们使用返回值代表是否有折射出现。
3. 递归追踪反射和折射
前篇在处理反射时，我们假设光线总是来自形状外部。然而，这里反射和折射都可以发生在形状的内部。因此，我们首先需要用 SDF 判断光线的发射坐标 $\mathbf{o}$ 位于形状之内（ $\phi(\mathbf{o}) < 0$ ）还是之外（ $\phi(\mathbf{o}) > 0$ ）。由于 SDF 会为正或负，在检测光线步进是否到形状表面时，我们要用绝对值来比较 $|\phi(\mathbf{o}+t \mathbf{d})|<\epsilon$ 。
我们为Result 结构加入形状的折射率：
typedef struct { float sd, emissive, reflectivity, eta; } Result;
如果光线从外至内，调用 refract() 时，传入 $1 /\eta$ ；从内至外则传入 $\eta$ 。
另外要注意的是，法线的方向需要是入射光的那一方，换句话说，当光线从形状内往外发射时，要反转法线方向。
float trace(float ox, float oy, float dx, float dy, int depth) {
float t = 1e-3f;
float sign = scene(ox, oy).sd > 0.0f ? 1.0f : -1.0f; // 内还是外？
for (int i = 0; i < MAX_STEP && t < MAX_DISTANCE; i++) {
float x = ox + dx * t, y = oy + dy * t;
Result r = scene(x, y);
if (r.sd * sign < EPSILON) { // 判断是否到达表面时要考虑内外
float sum = r.emissive;
if (depth < MAX_DEPTH && (r.reflectivity > 0.0f || r.eta > 0.0f)) {
float nx, ny, rx, ry, refl = r.reflectivity;
gradient(x, y, &nx, &ny);
nx *= sign; ny *= sign; // 在内的话，要反转法线
if (r.eta > 0.0f) {
if (refract(dx, dy, nx, ny, sign < 0.0f ? r.eta : 1.0f / r.eta, &rx, &ry))
sum += (1.0f - refl) * trace(x - nx * BIAS, y - ny * BIAS, rx, ry, depth + 1);
else
refl = 1.0f; // 全内反射
}
if (refl > 0.0f) {
reflect(dx, dy, nx, ny, &rx, &ry);
sum += refl * trace(x + nx * BIAS, y + ny * BIAS, rx, ry, depth + 1);
}
}
return sum;
}
t += r.sd * sign;
}
return 0.0f;
}
一些测试场景（矩形、凸透镜、凹透镜、半圆透镜）：
4. 结语
这篇仅加入了斯涅尔定律的应用，额外 20 行代码，就能渲染出类似玻璃材质的效果。我们在开始时假设了折射和反射的比率是恒定的（除了内全反射），这并不太接近真实的情况。我们在下篇会作改善。
本文的代码位于 refraction.c 及 m2.c
用 C 语言画光（六）：菲涅耳方程
在上一篇我们实现了折射，也假设了反射和折射的比重是恒定的。但现实上，这个比重与材质和入射角相关。
题图是我 2006 年在九寨沟拍摄的相片。相片底部，视角较接近垂直，可以看到较清晰的水下景物，倒影较暗；而在远处，视角较与水面平行，则几乎只看见倒影，难以看到水下景物。
你现在也可以把一个马克杯装满水，放在屏幕前观察水面。如果眼睛接近水平面，就可以清楚地看到屏幕倒影；如果视点往上移，那么屏幕倒影变暗，更多看到杯里的颜色。
这个现象称为菲涅耳反射（Fresnel reflectance）。奥古斯丁·菲涅耳（Augustin Fresnel， 1788－1827）为法国物理学家，注意 Fresnel 为法语姓氏，s 是不发音的。
1. 菲涅耳方程
菲涅耳方程（Fresnel equation）描述了光线经过两个介质的界面时，反射和透射的光强比重。参考下图的符号：
当中 $R \in [0, 1]$ 为反射比，因能量守恒，透射比为 $T = 1 - R$ 。
对于电介质（dielectric）而言，菲涅耳方程为：
$$ \begin{aligned} R_s &= \left( \frac{\eta_1\cos\theta_i-\eta_2\cos\theta_t}{\eta_1\cos\theta_i+\eta_2\cos\theta_t} \right) ^ 2 \\ R_p &= \left( \frac{\eta_1\cos\theta_t-\eta_2\cos\theta_i}{\eta_1\cos\theta_t+\eta_2\cos\theta_i} \right) ^ 2 \\ \end{aligned} \tag{1} $$$R_s$ 和 $R_p$ 分别对应入射光的 s 偏振（senkrecht polarized）和 p 偏振（parallel polarized）所造成的反射比。图形学中通常考虑光是无偏振的（unpolarized），也就是两种偏振是等量的，所以可以取其平均值：
$$ R = \frac{R_s + R_p}{2} \tag{2} $$那么，给定入射角、透射角的余弦，以及两个介质的折射率，就可简单实现：
float fresnel(float cosi, float cost, float etai, float etat) {
float rs = (etat * cosi - etai * cost) / (etat * cosi + etai * cost);
float rp = (etai * cosi - etat * cost) / (etai * cosi + etat * cost);
return (rs * rs + rp * rp) * 0.5f;
}
透射角实际上是使用斯涅尔定律从入射角得出的，所以假设折射率是常数，菲涅耳反射比是入射角的函数。下图展示了光线从空气射向不同材质时的菲涅耳反射比（来自[1] P. 233，$R$ 为菲涅耳反射比）：
在趋向 $90^{\circ}$ 入射角时，不论什么材质反射比都趋向 1。
另外，导体（conductor）和半导体（semiconductor）的菲涅耳方程会更复杂一些，而且不同对波长的影响较大（上图中红色和紫色曲线分别对应 RGB 三种波长），我们稍后再讨论。
2. 修改实现
这个函数用于取代原来恒定的 $\textbf{reflectivity}$ ，我们原来的追踪代码是这样的：
if (depth < MAX_DEPTH && (r.reflectivity > 0.0f || r.eta > 0.0f)) {
/* ... */
if (r.eta > 0.0f) {
if (refract(dx, dy, nx, ny, sign < 0.0f ? r.eta : 1.0f / r.eta, &rx, &ry))
sum += (1.0f - refl) * trace(x - nx * BIAS, y - ny * BIAS, rx, ry, depth + 1);
else
refl = 1.0f; // Total internal reflection
}
/* ... */
}
当有折射发生时，我们利用点积计算入射角和透射角的余弦，然后用 $\textbf{fresnel()}$ 获取反射比。和调用 $\textbf{refract()}$ 时相似，我们要考虑光是从空气（$\eta=1$）进入物体，还是相反：
if (depth < MAX_DEPTH && (r.reflectivity > 0.0f || r.eta > 0.0f)) {
/* ... */
if (r.eta > 0.0f) {
if (refract(dx, dy, nx, ny, sign < 0.0f ? r.eta : 1.0f / r.eta, &rx, &ry)) {
float cosi = -(dx * nx + dy * ny);
float cost = -(rx * nx + ry * ny);
refl = sign < 0.0f ? fresnel(cosi, cost, r.eta, 1.0f) : fresnel(cosi, cost, 1.0f, r.eta);
sum += (1.0f - refl) * trace(x - nx * BIAS, y - ny * BIAS, rx, ry, depth + 1);
}
else
refl = 1.0f; // Total internal reflection
}
/* ... */
}
我们比较一下，上一篇使用恒定反射比和本篇使用菲涅耳方程的渲染结果：
比较明显的区别是，光线垂直射向凸透镜和凹透镜时，使用恒定反比射会产生很亮的反射光斑，而使用菲涅耳方程的反射光斑就会暗得多，大部分光能穿过了透镜。
3. Schlick 近似
鉴于导体的菲涅耳方程较复杂，Schlick1 提供了一个近似的函数：
$$ R(\theta_i) \approx R(0) + (1-R(0))(1-\cos\theta_i)^5 \tag{3} $$此函数只需要对材质提供光线垂直反射时的 $R(0)$ 值，就能近似地计算出不同入射角旳菲涅耳反射比。下图展示此近似函数（虚线）与菲涅耳方程（实线）的比较（自 [1] P. 234）：
对于电介质，我们可以使用 (2) 计算出 $R(0)$ ：
$$ R(0) = \left( \frac{\eta_i - \eta_t}{\eta_i + \eta_t} \right)^2 \tag{4} $$对于导体，我们就需要提供 $R(0)$ 的值（可能需要多个频率，如 RGB）。另一点要注意的事，从低折射率材质到高折射率材质时，要使用 $\theta_t$ ：
float schlick(float cosi, float cost, float etai, float etat) {
float r0 = (etai - etat) / (etai + etat);
r0 *= r0;
float a = 1.0f - (etai < etat ? cosi : cost);
float aa = a * a;
return r0 + (1.0f - r0) * aa * aa * a;
}
这个函数采用和 fresnel() 相同的接口，可以在例子中替换。但如果要支持导体，则要改变接口。
这个例子中用 Schlick 近似并没有肉眼能分辨的差异，就不显示渲染结果了。
4. 结语
我们依据菲涅耳方程，以入射角及材质特性推断出光线的反射比和透射比，这样更接近真实世界的情况。
Schlick 近似被广泛应用在实时的基于物理渲染（physically based rendering, PBR）中，除了运算简单快捷，以 $R(0)$ 作为参数更容易让美术理解──白光直射材质所反射的颜色。用物理量（如本文未提及的波阻抗和磁导率）作为参数的话就很不友好。
本文的代码位于 fresnel.c。
参考
[1] Akenine-Möller, Tomas, Eric Haines, and Naty Hoffman. Real-time rendering, Third Edition. CRC Press, 2008.
Schlick, Christophe. An Inexpensive BRDF Model for Physically-based Rendering. Computer Graphics Forum, Vol. 13, No. 3, 1994. ↩︎