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

用 C 语言画光(五):折射

上篇反射(reflection)成对出现的光学现象就是折射(refraction)。当光从一个介质(medium)进入另一个介质,便会出现这两种现象。

基于能量守恒,当没有能量散失时,反射光的能量与折射光的能量之和等于入射光的能量。在本篇中我们先假设反射和折射的比率是恒定的,例如一个形状的表面总是会反射 20% 的入射光,并折射 80% 的入射光。

但我们知道,折射的方向会改变(所以才称为「折」射),那么怎样计算折射方向呢?

1. 斯涅尔定律

读者可能在中学阶段已学会斯涅尔定律(Snell's law),它描述了入射角、折射角与介质折射率(refractive index)的关系:

\begin{align} \eta_1 \sin \theta_1 &= \eta_2 \sin \theta_2\\ \Leftrightarrow \frac{\sin \theta_2}{\sin \theta_1} &= \frac{\eta_1}{\eta_2} \end{align}\\

\eta_1 > \eta_2 时,例如光线从玻璃( \eta_1\approx1.5 )进入空气( \eta_2\approx 1 ),会有机会出现内全反射(total internal reflection),从斯涅尔定律分析:

 \sin \theta_2 = \frac{\eta_1}{ \eta_2} \sin \theta_1\\

如果 \eta_1 > \eta_2 \Leftrightarrow \frac{\eta_1}{\eta_2} > 1 ,某些入射角 \theta_1 范围会令右侧大于 1,左侧便会超出正弦的值域 [-1, 1] 。在这种情况下,不会发生折射,所有入射光会完全反射。我们把上式左侧设为 1,则可求出临界角(critical angle):

1 = \frac{\eta_1}{ \eta_2} \sin \theta_c\\ \Leftrightarrow \theta_c = \sin^{-1}\frac{\eta_2}{ \eta_1}\\

2. 折射矢量

我们以矢量表示光线方向,需要求出折射方向。这会比上一篇反射复杂一点点。先看图:

与反射不一样,折射最好使用单位矢量表示入射方向 \hat{\mathbf{i}} 和折射方向 \hat{\mathbf{t}} 。这里还定义了界面的切线方向 \hat{\mathbf{m}} 。从上图可以看到,折射矢量可从法线方向与切线方向的矢量之和得出:

\hat{\mathbf{t}} = \sin \theta_2 \hat{\mathbf{m}} -\cos \theta_2 \hat{\mathbf{n}}\tag{1}

接下来,我们要把上式各部分以 \hat{\mathbf{i}}, \hat{\mathbf{n}},\eta_1,\eta_2 表示。切线 \hat{\mathbf{m}} 可以从入射方向 \hat{\mathbf{i}} 分解得出:

\begin{align} \hat{\mathbf{i}} &= (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})\hat{\mathbf{n}}\ + \sin \theta_1 \hat{\mathbf{m}}\\ \Leftrightarrow \hat{\mathbf{m}} &= \frac{\hat{\mathbf{i}} - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})\hat{\mathbf{n}}}{\sin \theta_1 }\tag{2} \end{align}

然后 \sin\theta_2 可通过斯涅尔定律以 \sin\theta_1, \eta_1, \eta_2 表示,而 \sin\theta_1 又可以转换成点积:

\begin{align} \sin \theta_2 &= \frac{\eta_1}{\eta_2} \sin \theta_1 = \frac{\eta_1}{\eta_2} \sqrt{1 - \cos^2\theta_1} = \frac{\eta_1}{\eta_2} \sqrt{1 - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})^2}\\ \end{align}\tag{3}

有正弦也就可以求余弦:

\begin{align} \cos\theta_2 &= \sqrt{1 - \sin^2\theta_2} = \sqrt{1 - \left(\frac{\eta_1}{\eta_2}\right)^2 \left(1 - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})^2\right)}\\ \end{align}\tag{4}

把(2)、(3)、(4) 代入 (1),简化:

\begin{align} \hat{\mathbf{t}} &= \sin \theta_2 \hat{\mathbf{m}} -\cos \theta_2 \hat{\mathbf{n}}\\ &= \frac{\sin\theta_2}{\sin\theta_1} \left( \hat{\mathbf{i}} - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})\hat{\mathbf{n}}\right) - \sqrt{1 - \left(\frac{\eta_1}{\eta_2}\right)^2 \left(1 - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})^2\right)} \hat{\mathbf{n}}\\ &= \frac{\eta_1}{\eta_2} \left( \hat{\mathbf{i}} - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})\hat{\mathbf{n}}\right) - \sqrt{1 - \left(\frac{\eta_1}{\eta_2}\right)^2 \left(1 - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})^2\right)} \hat{\mathbf{n}}\\ &= \frac{\eta_1}{\eta_2} \hat{\mathbf{i}} - \left(\frac{\eta_1}{\eta_2}(\hat{\mathbf{i}}\cdot\hat{\mathbf{n}}) + \sqrt{1 - \left(\frac{\eta_1}{\eta_2}\right)^2 \left(1 - (\hat{\mathbf{i}}\cdot\hat{\mathbf{n}})^2\right)} \right)\hat{\mathbf{n}}\\ \end{align}\\

噢,计算折射矢量竟然一个三角函数都不需要!由于上式只需折射率之比值,我们的函数使用 \eta =\eta_1 / \eta_2 作为参数,实现如下:

int refract(float ix, float iy, float nx, float ny, float eta, float* rx, float* ry) {
    float idotn = ix * nx + iy * ny;
    float k = 1.0f - eta * eta * (1.0f - idotn * idotn);
    if (k < 0.0f)
        return 0; // 全内反射
    float a = eta * idotn + sqrtf(k);
    *rx = eta * ix - a * nx;
    *ry = eta * iy - a * ny;
    return 1;
}

当根号内的值为负数,也是代表全内反射,我们使用返回值代表是否有折射出现。

3. 递归追踪反射和折射

前篇在处理反射时,我们假设光线总是来自形状外部。然而,这里反射和折射都可以发生在形状的内部。因此,我们首先需要用 SDF 判断光线的发射坐标 \mathbf{o} 位于形状之内( \phi(\mathbf{o})<0 )还是之外( \phi(\mathbf{o})>0 )。由于 SDF 会为正或负,在检测光线步进是否到形状表面时,我们要用绝对值来比较 |\phi(\mathbf{o}+t\mathbf{d}) |<\epsilon

我们为 \texttt{Result} 结构加入形状的折射率:

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

如果光线从外至内,调用 \texttt{refract()} 时,传入 1 / \eta ;从内至外则传入 \eta

另外要注意的是,法线的方向需要是入射光的那一方,换句话说,当光线从形状内往外发射时,要反转法线方向。

float trace(float ox, float oy, float dx, float dy, int depth) {
    float t = 1e-3f;
    float sign = scene(ox, oy).sd > 0.0f ? 1.0f : -1.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 * sign < EPSILON) { // 判断是否到达表面时要考虑内外
            float sum = r.emissive;
            if (depth < MAX_DEPTH && (r.reflectivity > 0.0f || r.eta > 0.0f)) {
                float nx, ny, rx, ry, refl = r.reflectivity;
                gradient(x, y, &nx, &ny);
                nx *= sign; ny *= sign; // 在内的话,要反转法线
                if (r.eta > 0.0f) {
                    if (refract(dx, dy, nx, ny, sign < 0.0f ? r.eta : 1.0f / r.eta, &rx, &ry))
                        sum += (1.0f - refl) * trace(x - nx * BIAS, y - ny * BIAS, rx, ry, depth + 1);
                    else
                        refl = 1.0f; // 全内反射
                }
                if (refl > 0.0f) {
                    reflect(dx, dy, nx, ny, &rx, &ry);
                    sum += refl * trace(x + nx * BIAS, y + ny * BIAS, rx, ry, depth + 1);
                }
            }
            return sum;
        }
        t += r.sd * sign;
    }
    return 0.0f;
}

一些测试场景(矩形、凸透镜、凹透镜、半圆透镜):

4. 结语

这篇仅加入了斯涅尔定律的应用,额外 20 行代码,就能渲染出类似玻璃材质的效果。我们在开始时假设了折射和反射的比率是恒定的(除了内全反射),这并不太接近真实的情况。我们在下篇会作改善。

本文的代码位于 refraction.cm2.c

编辑于 2017-12-11

文章被以下专栏收录