一步一步实现极简光线追踪

关于光线追踪实现的教程网上有很多，但是绝大部分教程都非常难理解，或者涉及很多第三方库，或者对其中的一些基础知识介绍很少，本文将采用纯软件来实现光线追踪算法，针对光线追踪的全过程进行详细的介绍，为了便于理解，抛弃窗口显示功能，将最终的结果输出到图片文件中。

程序的实现是基于 Unity3D 的，如果需要 C/C++ 版本，可以参考原教程，对于实现过程中涉及的一些几何计算，如果实在不理解或者不感兴趣，只需要清楚它的作用即可，不是说它不重要，相反，它是必不可少的，但是本文的重点是理解光线追踪算法，因此会尽可能的避免阻碍理解的内容。

0x01 输出图片

第一步是实现图片的存储输出功能，在实时渲染中，会将渲染结果显示到窗口视图，那样会做很多无关的工作，这里只是将渲染结果保存到图片中，减轻我们的理解压力。

在 Unity 中实现该功能非常的简单，实现代码在这里。

输出结果如下：

显示图形

只有渐变背景有点单调，下面在背景上画一个圆球，最终的效果如下：

0x02 定义球体

首先需要定义一个球体，代码在这里：

关于上面射线检测证明过程，可以参考这里，过程比较详细。

0x03 定义摄像机

定义好图形之后，如何将图形显示到我们的图片中呢？换句话说，如何确定图片的哪些像素显示球呢？为了实现这一功能，需要引入一个摄像机，想象在一个三维场景中，在一架摄像机前面放置这一个球，最终摄像机会记录下特定范围内的画面，如下图所示：

我们知道摄像机有一个拍摄范围，在计算机图形学中称这个范围为视锥体，为了表示这个范围，通常采用平截四棱锥来表示，具体的样子如下：

其中，靠近摄像机的截面叫做近裁剪面、远离摄像机的截面叫做远裁剪面。为了确定输出图片中每个像素值，我们可以把图片贴到近裁剪面，然后以相机为起点（这里放到原点位置）、向每个像素发射射线，检测与球是否相交，如果相交说明就需要显示，反之则不需要。

计算近裁剪面

那么下面问题是如何计算每个像素对应的三维坐标，因为要将图片贴到近裁剪面，所以首先需要计算一下近裁剪面的尺寸，这里又多了个概念叫视场角，它是指视锥体两边的夹角，根据这个角度可以计算出近裁剪面的高度，而宽度可以根据指定的宽高比$aspect$计算，如下图所示：

/// <summary>
/// 获取近裁剪面的尺寸
/// </summary>
/// <param name="near">近裁剪面距离摄像机的距离</param>
/// <param name="fov">视场角</param>
/// <param name="aspect">宽高比</param>
/// <returns></returns>
private Vector2 GetNearClipSize(float near, float fov, float aspect)
{
    float height = Mathf.Tan(fov * 0.5f) * 2.0f * near;
    float width = height * aspect;
    return new Vector2(width, height);
}

屏幕坐标转世界坐标

有了近裁剪面的尺寸之后，下一步就是怎么把宽度为width、高度为height的图片贴到近裁剪面（映射的一个过程），代码如下：

/// <summary>
/// 根据像素的索引，计算世界坐标
/// </summary>
/// <param name="i">像素的水平索引</param>
/// <param name="j">像素的垂直索引</param>
/// <returns></returns>
public Vector3 ScreenToWorld(int i, int j, int width, int height)
{
    float ix = i + 0.5f;    // 像素的水平中心
    float jy = j + 0.5f;    // 像素的垂直中心

    float aspect = (float)width / height;
    Vector2 nearClipSize = GetNearClipSize(1, 60.0f * Mathf.Deg2Rad, aspect);

    // 范围[0,nearClipSize.x]
    float x = ix * nearClipSize.x / width;
    float y = jy * nearClipSize.y / height;

    // 范围：[-nearClipSize.x / 2,nearClipSize.x / 2]
    x -= nearClipSize.x * 0.5f;
    y -= nearClipSize.y * 0.5f;

    return new Vector3(x, y, -1);
}

完整代码在这里。

0x04 渲染图形

经过一系列计算，终于能够根据像素获取到三位坐标了，修改步骤01中的 GetPixel 函数，实现在背景色(0.2,0.7,0.8)上显示一个球，实现步骤如下：

定义一个相机，并根据输入的像素位置，计算追踪射线；
定义一个球体；
进行射线检测，如果射线穿过球体，则返回球体的颜色，反之则返回背景色。

代码地址。

0x05 多添几个球

再向场景中添加几个不同颜色的球，最终的显示效果如下：

为了表示球的颜色信息，这里用材质来表示，同时定义 RenderObject 类来保存小球的使用的材质，查看代码。

当场景存在多个球时，就需要所有物体进行射线检测，为了方便计算，这里定义 World 对象来存储所有的球，同时为了方便查看各种效果，在场景中增加一个棋盘格平面，完整代码。

然后修改渲染部分代码：

创建几个材质；
创建 World 对象，并添加几个球体，为球体绑定材质；
修改 CastRay 方法，由单个球的射线检测，改为整个场景的检测。

完成代码在这里。

添加灯光

上面实现虽然我们绘制的是三维的球，但是结果看上去只是二维，这是因为没有光影的明暗变化，我们之所以能够看到各种颜色的物体，是因为光线照射到物体表面，经过反射后进入到人眼中，最终经过眼睛感光成像才看到物体。

在这里我们不用考虑人眼是如何感光成像的，只需要关心光线是经过怎样的过程进入到眼睛的，光线照射到物体的表面之后，会分为两部分，第一部分会在物体表面发生散射，第二部分会进入到物体内部，它们都会产生不同的视觉特性。

真实的光照是非常复杂的，因此在图形学中会采用经验模型来模拟光照的效果，这里采用标准光照模型来计算光照，它的基本方法是将进入人眼的光照分为四部分：

自发光
漫反射
高光反射
环境光

0x06 漫反射

当光照照射向物体表面时，因为物体的表面时凹凸不平的（这里的表面指的是微表面，可以认为是分子甚至更小的），所以会随机散射到各个方向，因为方向是随机，可以认为在任何方向上都是一样的，和人眼在什么位置是无关的，但是散射的强度与光线的入射角度相关。它们之间的关系符合兰伯特定律：反射光线的强度与表面法线和光源入射方向之间的夹角的余弦值成正比，计算公式如下：

$C_{diffuse} = C_{light} * m_{diffuse} * max(0,\vec{n}· \vec{l})$

最终效果如下：

先定义灯光，点光源的主要参数包括位置、光照强度，代码如下：

public class PointLight
{
    public Vector3 position;
    public float intensity;

    public PointLight(Vector3 pos, float i)
    {
        this.position = pos;
        this.intensity = i;
    }
}

修改渲染部分:

创建 PointLight 对象；
修改 CastRay 方法，当射线击中对象时，添加漫反射光照。

最终代码在这里。

0x07 高光反射

高光反射用于描述光线照射到物体表面时，在表面完全镜面反射方向的辐射量，这里采用 Phong 模型来计算高光，涉及的参数包括光照方向、表面法线、视线方向、镜面反射方向，如下图所示：

反射方向能够根据其他参数计算获得：

$\vec{r} = 2(\vec{n} · \vec{l})\vec{n} - \vec{l}$

根据 Phong 模型计算高光：

$C_{specular} = (C_{light} * m_{specular}) \max(0,\vec{v}· \vec{r})^{m_{gloss}}$

其中 $m{specular}$表示高光的颜色，$m{gloss}$表示反光度，该值越大则高光的亮点越小，修改代码如下：

Material 内增加高光属性；
修改 CastRay 函数，增加高光计算。

具体代码在这里。

0x08 添加阴影

为了生成阴影，就需要计算每个像素是否在阴影范围之内，因此就需要检测是否被其他物体遮住光线，检测射线的起始位置从检测点向法线反向偏移，方向指向光源位置，这样的射线检测到交叉物体之后，说明该检测位置被其他物体遮挡住了，此时只需跳过光照处理，起始位置计算方式如下：

1	orign = Vector3.Dot(ligthDir,normal) < 0? point - normal0.001f : point + normal0.001f;

实现代码在这里。

0x09 添加反射

前面提到，光线照射到物体表面时，会分为两部分，第一部分是进入到物体的内部，第二部分是反射到环境中，而这部分反射光照，同样会再次照射到其他物体之上，但是反射并不会无限进行下去，随着反射会不断的衰减，所以这里会设置一个深度值来表示反射的次数，反射位置的计算方式和阴影起始点计算方式相同，修改 CastRay 代码如下：

private Color CastRay(Vector3 origin, Vector3 dir, RenderObject[] objs, Light[] lights,int depth)
{
    Material material = null;
    Vector3 normal = Vector3.zero;
    Vector3 hitPoint = Vector3.zero;

    if (depth > 4 || !world.SceneIntersect(origin, dir, ref hitPoint, ref normal, ref material))
        return new Color(0.2f, 0.7f, 0.8f);
    else
    {
        Vector3 reflDir = -Reflect(dir, normal).normalized;
        Vector3 refl_Origin = Vector3.Dot(reflDir, normal) < 0 ? hitPoint - normal * 0.001f : hitPoint + normal * 0.001f;
        Color refl_Color = CastRay(refl_Origin, reflDir, depth);

        float diffuse_light_intensity = 0.0f;
        float specular_light_intensity = 0.0f;
        for (int i = 0; i < lights.Length; i++)
        {
            Vector3 lightDir = (lights[i].position - hitPoint).normalized;

            Vector3 sdRayOrigin = Vector3.Dot(normal, lightDir) < 0 ? hitPoint - normal * 0.001f : hitPoint + normal * 0.001f;
            Vector3 sdHitNormal = Vector3.zero;
            Vector3 sdHitPoint = Vector3.zero;
            Material sdHitMaterial = null;
            if (world.SceneIntersect(sdRayOrigin, lightDir, ref sdHitPoint, ref sdHitNormal, ref sdHitMaterial) &&
                Vector3.Distance(sdRayOrigin, sdHitPoint) < Vector3.Distance(hitPoint, lights[i].position))
            {// 检测在光源和物体之间，是否存在其他遮挡物
                continue;
            }

            diffuse_light_intensity += lights[i].intensity * Mathf.Max(0.0f, Vector3.Dot(lightDir, normal));
            specular_light_intensity += Mathf.Pow(Mathf.Max(0, Vector3.Dot(dir, Reflect(-lightDir, normal))), material.gloss) * lights[i].intensity;
        }

        return material.diffuse * diffuse_light_intensity +
                        material.specular * specular_light_intensity + 
                        refl_Color * material.reflectAmount;
    }
}

完整代码在这里。

0x0A 添加折射

无论漫反射、高光反射、镜面反射都属于光的反射部分，而进入物体的部分光线，会产生折射效果，其入射角和折射角之间可以用斯涅耳定律来表示。

其中 $\theta{1}$ 表示入射角，$\theta{2}$ 表示折射角，$n_1$ 是上层介质的折射率，$n_2$ 表示下层介质的折射率，$\vec{QO}$表示入射向量，$\vec{OP}$表示折射后光线的方向，$\vec{ON}$ 表示法线向量，以上这些向量都为单位向量。

求解折射向量 $\vec{OP}$ 过程如下：

分解向量$\vec{OP}$：

$\begin{aligned} \vec{OP} &= \vec{OK} + \vec{OB} \\ &= \vec{OM} * \sin{\theta_2} - \vec{ON} * \cos{\theta_2} \\ \end{aligned} \tag{1}$

分解 $\vec{QO}$ 向量：

$\begin{aligned} \vec{QO} &= \vec{OM} * \sin{\theta_1} + \vec{ON} * \cos{\theta_1} \\ \end{aligned} \tag{2}$

根据点乘 $\vec{QO}$、$\vec{ON}$:

$\begin{aligned} \vec{QO} \cdot \vec{ON} &= |\vec{QO}|\cdot |\vec{ON}| \cos{\theta_1} = \cos{\theta_1} \\ \sin{\theta_1} &= \sqrt{1 - {\cos^2{\theta_1}}} \\ \end{aligned} \tag{3}$

将 (3) 带入 (2) 式可得:

$\begin{aligned} \vec{OM} &= \frac{\vec{QO} - \vec{ON} * \cos{\theta_1}}{\sin{\theta_1}} \\ &= \frac{\vec{QO} - \vec{ON} * (\vec{QO} \cdot \vec{ON})}{\sin{\theta_1}} \\ \end{aligned} \tag{4}$

根据斯涅耳方程:

$\begin{aligned} n_1*\sin{\theta_1} &= n_2 * \sin{\theta_2} \\ \cos{\theta_2} &= \sqrt{1 - (\frac{n_1}{n_2})^2 * ({1-{\cos^2{\theta_1}}})} \\ &= \sqrt{1 - (\frac{n_1}{n_2})^2 * ({1-({\vec{QO}\cdot \vec{ON}})^2})} \\ \sin{\theta_2} &= \sqrt{1 - \cos^2{\theta_2}} \\ \end{aligned} \tag{5}$

将 (4)、(5) 带入 (1) 式可得：

$\begin{aligned} \vec{OP} &= \vec{OM} * \sin{\theta_2} - \vec{ON} * \cos{\theta_2} \\ &= (\frac{\vec{QO} - \vec{ON} * (\vec{QO} \cdot \vec{ON})}{\sin{\theta_1}} ) * \sin{\theta_2} - \vec{ON} * (\sqrt{1 - (\frac{n_1}{n_2})^2 * ({1-({\vec{QO}\cdot \vec{ON}})^2})}) \\ &= (\vec{QO} - \vec{ON} * (\vec{QO} \cdot \vec{ON})) * (\frac{n_1}{n_2}) - \vec{ON} * (\sqrt{1 - (\frac{n_1}{n_2})^2 * ({1-({\vec{QO}\cdot \vec{ON}})^2})}) \\ &= \vec{QO}*\frac{n_1}{n_2} - \vec{ON} *((\vec{QO}\cdot \vec{ON}) * \frac{n_1}{n_2} + (\sqrt{1 - (\frac{n_1}{n_2})^2 * ({1-({\vec{QO}\cdot \vec{ON}})^2})})) \\ \end{aligned}$

上式中需要开平方根，所以需要其下的值不小于0才有解，当其下值小于0时，上式无解，即不存在折射向量。

有了上面的计算结果之后，代码实现就非常简单了，具体代码如下：

private Vector3 Refract(Vector3 Q, Vector3 N, float eta_t, float eta_i = 1.0f)
{
  float cos_theta1 = Vector3.Dot(-Q, N);
  if (cos_theta1 < 0)
      return Refract(Q, -N, eta_i, eta_t);

  float eta = eta_i / eta_t;
  float k = 1 - eta * eta * (1 - cos_theta1 * cos_theta1);
  return k < 0 ? new Vector3(0, 0, 0) : (Q * eta - N * (eta * cos_theta1 + Mathf.Sqrt(k)));
}

完整代码在这里。

总结

再次把射线追踪的示意图放出来：

整个过程是从眼睛（Camera）开始，向三维场景中发出射线，获取碰撞点的颜色，而这个颜色可能有一下几种情况：

没有碰到任何物体，返回背景色；
碰到物体，如果存在灯光，则需要根据光照模型进行计算，结果由几部分组成：
1. 漫反射：与光线、法线有关；
2. 高光反射：与光线、视线、法线有关；
3. 阴影：从碰撞点向光源发出射线，检测是否存在遮挡；
4. 反射：反射的哪个物体？需要计算射线的反射向量，再次进行射线检测；
5. 折射：折射的哪个物体？需要计算射线的折射向量，再次进行射线检测。