用 C 语言画光

野性酷女 2022-06-06 12:23 65阅读 0赞

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

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

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

## 1. 光 ##

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

![F(x, y) = \\int\_0^\{2\\pi\} L(x, y, \\theta)\\mathrm\{d\}\\theta][F_x_ y_ _ _int_0_2_pi_ L_x_ y_ _theta_mathrm_d_theta]

当中 ![L(x, y, \\theta)][L_x_ y_ _theta] 代表在二维坐标 ![(x, y)][x_ y] 在 ![\\theta][theta] 方向有多少光经过。

我们想采样这个 ![F(x, y)][F_x_ y] 存储在图像中，并利用 [svpng()][svpng] 输出结果，那么我们的 ![\\texttt\{main()\}][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\]][x_ y_ _in _0_ 1_ _times _0_ 1] 。我们把结果映射至 ![\\\{0, 1, \\ldots, 255\\\}][0_ 1_ _ldots_ 255] 。

## 2. 蒙地卡罗积分 ##

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

在这个问题中，我们随机采样 ![N][] 个方向 ![\{\\theta\_1, \\theta\_2, \\ldots, \\theta\_N\}][theta_1_ _theta_2_ _ldots_ _theta_N] ，然后计算 ![L(x, y, \\theta\_i)][L_x_ y_ _theta_i] 的平均值：

![F(x, y) \\approx \\frac\{2\\pi\}\{N\}\\sum\_\{i=1\}^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)\}][texttt_trace_ox_ oy_ dx_ dy] 函数代表从 ![\\mathbf\{o\}][mathbf_o] 位置从单位矢量 ![\\hat\{\\mathbf\{d\}\}][hat_mathbf_d] 方向接收到的光。

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

## 3. 光线步进 ##

通常，我们可以用光线追踪（ray tracing）方法，求出光线 ![\\mathbf\{r\}(t) = \\mathbf\{o\} + \\hat\{\\mathbf\{d\}\}t][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 _ _mathbb_R_2 _rightarrow _mathbb_R] 表示

1.  当 ![\\phi(\\mathbf\{x\}) > 0][phi_mathbf_x_ _ 0] ，表示坐标 ![\\mathbf\{x\}][mathbf_x] 位于场景形状之外，且 ![\\mathbf\{x\}][mathbf_x] 与最近形状边界的距离为 ![\\phi(\\mathbf\{x\})][phi_mathbf_x] ；
2.  当 ![\\phi(\\mathbf\{x\}) < 0][phi_mathbf_x_ _ 0 1] ，表示坐标 ![\\mathbf\{x\}][mathbf_x] 位于场景形状之内，且 ![\\mathbf\{x\}][mathbf_x] 与最近形状边界的距离为 ![\-\\phi(\\mathbf\{x\})][-_phi_mathbf_x] ；
3.  当 ![\\phi(\\mathbf\{x\})= 0][phi_mathbf_x_ 0] ，说明 ![\\mathbf\{x\}][mathbf_x] 刚好在形状边界上。

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

![\\phi\_\\text\{circle\}(\\mathbf\{x\}) = \\left \\| \\mathbf\{x\} - \\mathbf\{c\} \\right \\|-r][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 _ 0] ，逐步增加 ![t][] ，当 ![\\phi(\\mathbf\{r\}(t)) \\le 0][phi_mathbf_r_t_ _le 0] 代表我们到达到场景中某个形状的表面或内部。那么我们每次可以步进多远？由于 ![\\phi(\\mathbf\{x\})>0][phi_mathbf_x_0] 时，代表 ![\\mathbf\{x\}][mathbf_x] 距最近形状的距离，所以我们至少可以步进该距离而不会碰到任何形状！

![v2-17fd23c7b0af869f90bc2d374e56aa1e_hd.jpg][]图源 https://developer.nvidia.com/gpugems/GPUGems2/gpugems2\_chapter08.html

设场景只有一个发光的圆形，圆心为 ![(0.5, 0.5)][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\}][epsilon _ 10_-6] 表示足够近的阈值。整个程序完成，运算结果如下：

![v2-0a31e467fd69a2f08e0fcea6a1979ff2_hd.jpg][]

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

![v2-8d86325bc5b25ae3f9ac26b60547410a_hd.jpg][]

然而，随着 ![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;             // 分层采样
            // ...

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

![v2-44cbaf8810003bb099f6894c0c44810a_hd.jpg][]

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

我们可以结合上面两种采样方式，就是先分为 ![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; // 抖动采样

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

![v2-6d1dba9dd3541e0c1b17d2211923a559_hd.jpg][]

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

![v2-345085c75f4e84db7dde6bdf383f4d47_hd.jpg][]

## 5. 结语 ##

这个完整程序位于 [basic.c][]。如果不含加载头文件及常数定义，仅有 30 行代码。你可以改一下圆形位置、大小、光度，测试时可用较小的 W、H 和 N 加速运行过程。

本文简单介绍了蒙地卡罗积分、光线步进、分层采样等概念。下一篇《[构造实体几何][Link 5]》会讲如何在场景中加入更多形状，将会显示出阴影。而之后也会尝试实现一些光学法则，生成更有趣的图形。

## 参考 ##

\[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

[Link 1]: https://www.zhihu.com/question/20187195/answer/34873279
[Link 2]: https://www.zhihu.com/question/27015321/answer/35028446
[Link 3]: https://zhuanlan.zhihu.com/p/24688522
[Link 4]: https://zhuanlan.zhihu.com/p/30553006
[miloyip_light2d]: http://link.zhihu.com/?target=https%3A//github.com/miloyip/light2d
[x_ y]: http://www.zhihu.com/equation?tex=%28x%2C+y%29
[F_x_ y_ _ _int_0_2_pi_ L_x_ y_ _theta_mathrm_d_theta]: http://www.zhihu.com/equation?tex=F%28x%2C+y%29+%3D+%5Cint_0%5E%7B2%5Cpi%7D+L%28x%2C+y%2C+%5Ctheta%29%5Cmathrm%7Bd%7D%5Ctheta
[L_x_ y_ _theta]: http://www.zhihu.com/equation?tex=L%28x%2C+y%2C+%5Ctheta%29
[theta]: /images/20220606/b57793d50c454a6ea9dfc4504119b842.png
[F_x_ y]: http://www.zhihu.com/equation?tex=F%28x%2C+y%29
[svpng]: https://zhuanlan.zhihu.com/p/26525083
[texttt_main]: /images/20220606/3f3539f29b794fc08458be17a6a1c9e9.png
[x_ y_ _in _0_ 1_ _times _0_ 1]: http://www.zhihu.com/equation?tex=%28x%2C+y%29+%5Cin+%5B0%2C+1%5D+%5Ctimes+%5B0%2C+1%5D
[0_ 1_ _ldots_ 255]: http://www.zhihu.com/equation?tex=%5C%7B0%2C+1%2C+%5Cldots%2C+255%5C%7D
[N]: /images/20220606/1d4af58263044b6291ef8680d225ee56.png
[theta_1_ _theta_2_ _ldots_ _theta_N]: http://www.zhihu.com/equation?tex=%7B%5Ctheta_1%2C+%5Ctheta_2%2C+%5Cldots%2C+%5Ctheta_N%7D
[L_x_ y_ _theta_i]: http://www.zhihu.com/equation?tex=L%28x%2C+y%2C+%5Ctheta_i%29
[F_x_ y_ _approx _frac_2_pi_N_sum_i_1_N L_x_ y_ _theta_i]: http://www.zhihu.com/equation?tex=F%28x%2C+y%29+%5Capprox+%5Cfrac%7B2%5Cpi%7D%7BN%7D%5Csum_%7Bi%3D1%7D%5EN+L%28x%2C+y%2C+%5Ctheta_i%29
[texttt_trace_ox_ oy_ dx_ dy]: http://www.zhihu.com/equation?tex=%5Ctexttt%7Btrace%28ox%2C+oy%2C+dx%2C+dy%29%7D
[mathbf_o]: http://www.zhihu.com/equation?tex=%5Cmathbf%7Bo%7D
[hat_mathbf_d]: /images/20220606/150d22c541e2466b941ea8ce6799bf6b.png
[2_pi]: /images/20220606/670c93b101584764934497d4ea4228bd.png%5Cpi
[Bimos]: /images/20220606/28b630f6297a46e481db1b07d4e3416b.png
[mathbf_r_t_ _ _mathbf_o_ _ _hat_mathbf_d_t]: http://www.zhihu.com/equation?tex=%5Cmathbf%7Br%7D%28t%29+%3D+%5Cmathbf%7Bo%7D+%2B+%5Chat%7B%5Cmathbf%7Bd%7D%7Dt
[phi _ _mathbb_R_2 _rightarrow _mathbb_R]: http://www.zhihu.com/equation?tex=%5Cphi+%3A+%5Cmathbb%7BR%7D%5E2+%5Crightarrow+%5Cmathbb%7BR%7D
[phi_mathbf_x_ _ 0]: /images/20220606/cc409308b5d546c582e9d73d2b18fa5c.png+%3E+0
[mathbf_x]: /images/20220606/d343c9308f054908893c0438dbb7981b.png
[phi_mathbf_x]: /images/20220606/cc409308b5d546c582e9d73d2b18fa5c.png
[phi_mathbf_x_ _ 0 1]: /images/20220606/cc409308b5d546c582e9d73d2b18fa5c.png+%3C+0
[-_phi_mathbf_x]: /images/20220606/951471fede284925a864588b140a2132.png
[phi_mathbf_x_ 0]: /images/20220606/cc409308b5d546c582e9d73d2b18fa5c.png%3D+0
[mathbf_c]: http://www.zhihu.com/equation?tex=%5Cmathbf%7Bc%7D
[r]: /images/20220606/8f1598c781c949bbb113c0d7f167af42.png
[phi_text_circle_mathbf_x_ _ _left _ _mathbf_x_ - _mathbf_c_ _right _-r]: http://www.zhihu.com/equation?tex=%5Cphi_%5Ctext%7Bcircle%7D%28%5Cmathbf%7Bx%7D%29+%3D+%5Cleft+%5C%7C+%5Cmathbf%7Bx%7D+-+%5Cmathbf%7Bc%7D+%5Cright+%5C%7C-r
[t _ 0]: http://www.zhihu.com/equation?tex=t+%3D+0
[t]: http://www.zhihu.com/equation?tex=t
[phi_mathbf_r_t_ _le 0]: http://www.zhihu.com/equation?tex=%5Cphi%28%5Cmathbf%7Br%7D%28t%29%29+%5Cle+0
[phi_mathbf_x_0]: /images/20220606/cc409308b5d546c582e9d73d2b18fa5c.png%3E0
[v2-17fd23c7b0af869f90bc2d374e56aa1e_hd.jpg]: /images/20220606/5629ad15b7b94563a13142cc218c0db4.png
[0.5_ 0.5]: http://www.zhihu.com/equation?tex=%280.5%2C+0.5%29
[0.1]: http://www.zhihu.com/equation?tex=0.1
[2]: /images/20220606/670c93b101584764934497d4ea4228bd.png
[epsilon _ 10_-6]: http://www.zhihu.com/equation?tex=%5Cepsilon+%3D+10%5E%7B-6%7D
[v2-0a31e467fd69a2f08e0fcea6a1979ff2_hd.jpg]: /images/20220606/980cc3b7f2b64215aa20eb2a35ef755f.png
[v2-8d86325bc5b25ae3f9ac26b60547410a_hd.jpg]: /images/20220606/d01b6fb3495842439b6e17bb8eef22fd.png
[v2-44cbaf8810003bb099f6894c0c44810a_hd.jpg]: /images/20220606/91b5b58e010546718dd4b6e355c6ee38.png
[v2-6d1dba9dd3541e0c1b17d2211923a559_hd.jpg]: /images/20220606/58879fcc80a8499b95436124aceb620d.png
[N_64]: /images/20220606/1d4af58263044b6291ef8680d225ee56.png%3D64
[v2-345085c75f4e84db7dde6bdf383f4d47_hd.jpg]: /images/20220606/87952b1106df4f69841344b1c1f8c254.png
[basic.c]: http://link.zhihu.com/?target=https%3A//github.com/miloyip/light2d/blob/master/basic.c
[Link 5]: https://zhuanlan.zhihu.com/p/30748318