用 C 语言画光

来源：互联网发布：linux设置防火墙端口编辑：程序博客网时间：2024/05/20 20:05

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

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

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

1. 光

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

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

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

我们想采样这个 $F(x, y)$ 存储在图像中，并利用 svpng() 输出结果，那么我们的 $\texttt{main()}$ 函数简单为：

#include "svpng.inc"#include <math.h>#define W 512#define H 512unsigned 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)$

代码实现也很浅白（设每像素向 64 个随机方向均匀采样／uniform sampling）：

#define TWO_PI 6.28318530718f#define N 64float 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}$ 表示

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

例如，圆心为 $\mathbf{c}$ 、半径为 $r$ 的圆形 SDF 定义为（更精确的说法是圆盘／disk）：

$\phi_\text{circle}(\mathbf{x}) = \left \| \mathbf{x} - \mathbf{c} \right \|-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(\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-6ffloat 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.

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

阅读全文

0 0