当前位置: 开发笔记 > 编程语言 > 正文

用C语言画光

作者：北极光的悲伤 | 来源：互联网 | 2023-10-11 17:38

之前写过用C语言画心形、圣诞树、雪花、直线，这次我们试玩光学。在这系列文章中，我们会在二维的画布上「绘画」一些基于物理的光。但作为趣味性的编程文章，我尽量用直觉、简单的方式去生成这些图像

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

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

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

1. 光

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

无论图像长宽是多少，这个二维空间的坐标范围都是 $(x, y) \in [0, 1] \times [0, 1]$ 。我们把结果映射至 $\{0, 1, \ldots, 255\}$ 。

2. 蒙地卡罗积分

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

在这个问题中，我们随机采样个方向 ${\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 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}$ 表示

当 $\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}$ 、半径为的圆形 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 ，逐步增加，当 $\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 ，它每点都向各方向发射个单位的光。那么光线步进可实现为：

#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}$ 表示足够近的阈值。整个程序完成，运算结果如下：

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

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

同样使用每像素 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