首发于Milo的编程
用 C 语言画光(四):反射

用 C 语言画光(四):反射

上一篇我们讨论了一些形状,本篇回到基于光学的模拟,先谈较简单的──光的反射(reflection)。我们假设形状的边界是平滑的,会产生镜面反射(specular reflection)。实现反射需要一点几何上的推导,以及考虑如何用 SDF 获得边界法线,最后把追踪过程改为递归方式进行。

1. 反射矢量

给定一个入射方向 \mathbf i 及边界法线 \hat {\mathbf n} ,我们可以根据下图计算出反射方向 \mathbf r

首先,把 \mathbf i 用点积投影至法线方向,得到法线方向的长度(注意法线需为单位矢量),然后再把这长度往法线方向缩放,得到图中蓝色虚线矢量 (\mathbf{i} \cdot \hat{\mathbf{n}})\hat{\mathbf{n}} 。之后,把 \mathbf i 减去两个蓝色虚线矢量就能得到反射矢量 \mathbf{r}

\mathbf r = \mathbf{i} - 2(\mathbf{i} \cdot \hat{\mathbf{n}})\hat{\mathbf{n}}\\

用 C 语言实现:

void reflect(float ix, float iy, float nx, float ny, float* rx, float* ry) {
    float idotn2 = (ix * nx + iy * ny) * 2.0f;
    *rx = ix - idotn2 * nx;
    *ry = iy - idotn2 * ny;
}

但问题是,我们怎样从场景中取得边界法线?

2. 梯度与法线

在此系列文章中,我们采用了 SDF 来构建场景。在形状的边界,距离场变化最大的方向便是法线方向。根据矢量微积分(vector calculus),一个纯量场(scalar field)的最大变化方向就是其梯度(gradient),所以我们需要求形状边界位置的 SDF 梯度:

\nabla \phi = \left ( \frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y} \right )\\

这梯度其实就是 SDF 分别在 xy 方向的变化率。SDF 几乎处处可导,而且,SDF 为程函方程的(Eikonal equation)特例,具特性:

\| \nabla \phi\|=1\\

换句话说,SDF 的梯度是单位矢量。

假设我们不知道 SDF 梯度的分析解,可利用数值方法,如中心差分法(central difference),求出梯度的近似值:

\begin{align} \nabla \phi &= \left( \frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y} \right)\\ &\approx \left( \frac{\phi(x + h, y) - \phi(x - h, y)}{2h}, \frac{\phi(x, y + h) - \phi(x, y - h)}{2h} \right) \end{align}\\

用 C 语言实现:

void gradient(float x, float y, float* nx, float* ny) {
    *nx = (scene(x + EPSILON, y).sd - scene(x - EPSILON, y).sd) * (0.5f / EPSILON);
    *ny = (scene(x, y + EPSILON).sd - scene(x, y - EPSILON).sd) * (0.5f / EPSILON);
}

我们用一个测试场景去看看它的梯度,当中有一个圆形、两个矩形:

Result scene(float x, float y) {
    Result a = {  circleSDF(x, y, 0.4f,  0.2f, 0.1f), 2.0f };
    Result b = {     boxSDF(x, y, 0.5f,  0.8f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f };
    Result c = {     boxSDF(x, y, 0.8f,  0.5f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f };
    return unionOp(unionOp(a, b), c);
}

我们改一下 \texttt{main()} 函数,用颜色把梯度可视化,红色通道显示 \frac{\partial \phi}{\partial x} ,绿色通道显示 \frac{\partial \phi}{\partial y}

for (int y = 0; y < H; y++)
    for (int x = 0; x < W; x++, p += 3) {
         float nx, ny;
         gradient((float)x / W, (float)y / H, &nx, &ny);
         p[0] = (int)((fmaxf(fminf(nx, 1.0f), -1.0f) * 0.5f + 0.5f) * 255.0f);
         p[1] = (int)((fmaxf(fminf(ny, 1.0f), -1.0f) * 0.5f + 0.5f) * 255.0f);
         p[2] = 0;
     }

上图中可以看到,大部分位置的梯度是连续的,至少我们现在需要在形状边界上的位置是连续的。

(注:由于有数值误差,计算出来的梯度分量可能略超出 [-1, 1] 的范围,所以在输出时限制其范围。但这个轻微的误差不影响我们计算反射。)

3. 递归追踪

我们考虑形状边界的反射率(reflectivity)为 [0, 1] 范围,为此增加一个字段:

typedef struct { float sd, emissive, reflectivity; } Result;

Result scene(float x, float y) {
    Result a = {  circleSDF(x, y, 0.4f,  0.2f, 0.1f), 2.0f, 0.0f };
    Result b = {     boxSDF(x, y, 0.5f,  0.8f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f, 0.9f };
    Result c = {     boxSDF(x, y, 0.8f,  0.5f, TWO_PI / 16.0f, 0.1f, 0.1f), 0.0f, 0.9f };
    return unionOp(unionOp(a, b), c);
}

如果反射率超过 1,总能量就变多,不符合能量守恒。少于 1 代表形状吸收了能量。

之后,修改 \texttt{trace()} ,当我们在向一个方向追踪光线时,若碰到会反射的形状时,就计算该边界相交点的法线,然后递归追踪反射方向,把收到的所有光加总。为免无限反射做成无限递归,需要限制反射次数:

#define BIAS 1e-4f
#define MAX_DEPTH 3

float trace(float ox, float oy, float dx, float dy, int depth) {
    float t = 0.0f;
    for (int i = 0; i < MAX_STEP && t < MAX_DISTANCE; i++) {
        float x = ox + dx * t, y = oy + dy * t;
        Result r = scene(x, y);
        if (r.sd < EPSILON) {
            float sum = r.emissive;
            if (depth < MAX_DEPTH && r.reflectivity > 0.0f) {
                float nx, ny, rx, ry;
                gradient(x, y, &nx, &ny);
                reflect(dx, dy, nx, ny, &rx, &ry);
                sum += r.reflectivity * trace(x + nx * BIAS, y + ny * BIAS, rx, ry, depth + 1);
            }
            return sum;
        }
        t += r.sd;
    }
    return 0.0f;
}

float sample(float x, float y) {
    float sum = 0.0f;
    for (int i = 0; i < N; i++) {
        float a = TWO_PI * (i + (float)rand() / RAND_MAX) / N;
        sum += trace(x, y, cosf(a), sinf(a), 0); // 开始时深度为 0
    }
    return sum / N;
}

另外,由于我们要向反射方向追踪,如果用原来的相交点,追踪时就为立即遇到 \texttt{r.sd < EPSILON} 而停止,所以我们稍微把相交点往法线方向偏移 \texttt{BIAS} 的距离,只要 \texttt{BIAS} > \texttt{EPSILON} ,就可以避免这个问题。但太大的话也会造成误差。结果:

我们可以看到两个矩形之间有第二次的反射。

利用平面和圆形的相对补集,就能做出凹面镜(concave mirror),弧形会产生题图中的焦散(caustics)现象:

Result scene(float x, float y) {
    Result a = { circleSDF(x, y, 0.4f, 0.2f, 0.1f), 2.0f, 0.0f };
    Result d = {  planeSDF(x, y, 0.0f, 0.5f, 0.0f, -1.0f), 0.0f, 0.9f };
    Result e = { circleSDF(x, y, 0.5f, 0.5f, 0.4f), 0.0f, 0.9f };
    return unionOp(a, subtractOp(d, e));
}

4. 结语

也许大家在中学物理中都学过镜面反射的知识,但用简单的程序实现出来总会令人觉得更有趣。

本文的重点是用数值方法求 SDF 在边界上的法线,以及通过递归来实现多次反射。

下一篇将讨论折射,那么,我们就可以画出穿过透镜的光了。

本文的源代码位于 reflection.c

编辑于 2017-11-19

文章被以下专栏收录