Milo的编程
首发于Milo的编程
用 C 语言画光(一):基础

用 C 语言画光(一):基础

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

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

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

1. 光

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

F(x, y) = \int_0^{2\pi} L(x, y, \theta)\mathrm{d}\theta\tag{1}

当中 L(x, y, \theta) 代表在二维坐标 (x, y)\theta 方向有多少光经过。

我们想采样这个 F(x, y) 存储在图像中,并利用 svpng() 输出结果,那么我们的 \texttt{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) 的平均值:

F(x, y) \approx \frac{2\pi}{N}\sum_{i=1}^N L(x, y, \theta_i)\tag{2}

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

当中 \texttt{trace(ox, oy, dx, dy)} 函数代表从 \mathbf{o} 位置从单位矢量 \hat{\mathbf{d}} 方向接收到的光。

(更新:本文不考虑实际单位,所以实现时把系数 2\pi 去掉了。感谢 @Bimos 提醒)

3. 光线步进

通常,我们可以用光线追踪(ray tracing)方法,求出光线 \mathbf{r}(t) = \mathbf{o} + \hat{\mathbf{d}}t 与场景的最近点。

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

为简单起见,本文采用光线步进(ray marching)方法(又称为球体追踪/sphere tracing [1]),场景只需要以带符号距离场(signed distance field, SDF) \phi : \mathbb{R}^2 \rightarrow \mathbb{R} 表示

  1. \phi(\mathbf{x}) > 0 ,表示坐标 \mathbf{x} 位于场景形状之外,且 \mathbf{x} 与最近形状边界的距离为 \phi(\mathbf{x})
  2. \phi(\mathbf{x}) < 0 ,表示坐标 \mathbf{x} 位于场景形状之内,且 \mathbf{x} 与最近形状边界的距离为 -\phi(\mathbf{x})
  3. \phi(\mathbf{x})= 0 ,说明 \mathbf{x} 刚好在形状边界上。

例如,圆心为 \mathbf{c} 、半径为 r 的圆形 SDF 定义为(更精确的说法是圆盘/disk):

\phi_\text{circle}(\mathbf{x}) = \left \| \mathbf{x} - \mathbf{c} \right \|-r\tag{3}

用代码表示:

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(\mathbf{r}(t)) \le 0 代表我们到达到场景中某个形状的表面或内部。那么我们每次可以步进多远?由于 \phi(\mathbf{x})>0 时,代表 \mathbf{x} 距最近形状的距离,所以我们至少可以步进该距离而不会碰到任何形状!

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

设场景只有一个发光的圆形,圆心为 (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;
}

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

可以看到图像中有很多噪点,这是由于蒙地卡罗积分法具随机性,计算出来的估值与精确数值会有误差。增加采样数目 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.

编辑于 2018-01-03

文章被以下专栏收录