用 C 语言画光

之前写过用 C 语言画心形、圣诞树、雪花、直线，这次我们试玩光学。

在这系列文章中，我们会在二维的画布上「绘画」一些基于物理的光。但作为趣味性的编程文章，我尽量用直觉、简单的方式去生成这些图像，仅介绍一些概念，和一般的正式计算机图形学内容不同。

本系列的源代码位于 miloyip/light2d。

1. 光

假设我们在一个二维的世界，这里有些会发光的二维形状，并暂时只考虑单色光。我们想知道的是，在这个空间中，每一点从 360 度各方向共有多少光经过。换成数学方式表示，我们想对每个二维坐标  求积分：

当中  代表在二维坐标  在  方向有多少光经过。

我们想采样这个  存储在图像中，并利用 svpng() 输出结果，那么我们的 函数简单为：

#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);
}

无论图像长宽是多少，这个二维空间的坐标范围都是  。我们把结果映射至  。

2. 蒙地卡罗积分

由于无法以解析式求解这个积分，我们使用蒙地卡罗积分法（Monte Carlo integration）。

在这个问题中，我们随机采样  个方向  ，然后计算  的平均值：

代码实现也很浅白（设每像素向 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;
}

当中  函数代表从  位置从单位矢量  方向接收到的光。

（更新：本文不考虑实际单位，所以实现时把系数  去掉了。感谢 @Bimos 提醒）

3. 光线步进

通常，我们可以用光线追踪（ray tracing）方法，求出光线  与场景的最近点。

然而，我们需要为每种几何形状编写与光线的相交函数（通常比较复杂），之后做一些效果可能还要提供相交点的法线（normal vector）。

为简单起见，本文采用光线步进（ray marching）方法（又称为球体追踪／sphere tracing [1]），场景只需要以带符号距离场（signed distance field, SDF）  表示

1. 当  ，表示坐标  位于场景形状之外，且  与最近形状边界的距离为  ；
2. 当  ，表示坐标  位于场景形状之内，且  与最近形状边界的距离为  ；
3. 当  ，说明  刚好在形状边界上。

例如，圆心为  、半径为  的圆形 SDF 定义为（更精确的说法是圆盘／disk）：

用代码表示：

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;
}

而所谓光线步进，就是我们从光线起始点  ，逐步增加  ，当  代表我们到达到场景中某个形状的表面或内部。那么我们每次可以步进多远？由于  时，代表  距最近形状的距离，所以我们至少可以步进该距离而不会碰到任何形状！

图源 https://developer.nvidia.com/gpugems/GPUGems2/gpugems2_chapter08.html

设场景只有一个发光的圆形，圆心为  ，半径为  ，它每点都向各方向发射 个单位的光。那么光线步进可实现为：

#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;
}

如果光线超过指定距离，或是步数太多，都终止步进。注意我们只能尽量接近表面，所以用  表示足够近的阈值。整个程序完成，运算结果如下：

可以看到图像中有很多噪点，这是由于蒙地卡罗积分法具随机性，计算出来的估值与精确数值会有误差。增加采样数目  ，就能令结果更精确（准确地说是减少方差／variance）：

然而，随着  上升，运算时间也线性上升。有没有方法可以改善？

4. 分层、抖动采样

形成噪点的原因在于，每个像素所得到的随机方向都很不一样。那么，如果我们不用随机方向，而是平分 360 度向  个方向采样，效果会如何？这种方式称为分层采样（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;             // 分层采样
        // ...

改变这一行代码的结果是：

很好，没有噪点，但也太过规律了！

我们可以结合上面两种采样方式，就是先分为  个角度区域，再从区域中均匀采样，这称为抖动采样（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; // 抖动采样

改变这一行代码的结果是：

同样使用每像素  次采样，仅仅改变采样方式（一行代码），效果就好很多：

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.

https://zhuanlan.zhihu.com/p/30745861