【笔记】不想玩海龟画图？用 C++ 来画画！

推荐自己跟着写自己的代码食用，而不是复制粘贴。=)

这篇文章基本来自《Ray Tracing in One Weekend》。

你的第一幅作品

ppm 图片格式

ppm 是一个简单的图像格式：

ppm 格式是一种图片格式，全名 Portable PixelMap。PPM格式的图片如果用文本编辑器打开会发现都是 ASCII 码而不是二进制。所以这种格式的图片方便操作，但是文件大小可能不是很友好。

ppm 的格式如下：

P3
m n
255
r g b
r g b
...

其中 $m,n$ 分别为图片宽度和高度，随后的 $r,g,b$ 为从上往下，从左往右，每个像素的颜色。

Hello World!

那么就尝试写一个代码试试吧：

#include <iostream>

int main() {
    int n = 256, m = 256;
    std::cout << "P3\n" << m << " " << n << "\n255\n";
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j) {
            int r = 0;
            int g = i;
            int b = j;
            std::cout << r << " " << g << " " << b << "\n";
        }
    return 0;
}

跑一下，把输出保存成 xxx.ppm，然后打开，应该会长这样：

如果你无法打开 ppm 格式的图片文件，可以在 VSCode 装 simple ppm viewer 插件，或者在线查看器，或者使用如下 python 程序（python <python-file-name> <input-image> <output-image>）把 ppm 格式的图片文件导出成 jpg 或 png（deepseek 写的）：

from PIL import Image
import argparse

def convert_ppm(input_path, output_path):
    try:
        img = Image.open(input_path)
        img.save(output_path)
    except Exception as e:
        raise RuntimeError(f"Failed to transform: {e}")

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("input")
    parser.add_argument("output")
    args = parser.parse_args()

    convert_ppm(args.input, args.output)

简单图像渲染

为了方便，请你掏出你尘封已久的三维计算几何的三维向量类。这里假设你已经实现好了三维向量类并且有一些基本数学计算的功能。

另外，我们还需要一个 color 类，为了方便偷懒，我们钦定颜色类的三个参数取值为 $[0,1]$ 表示红色、蓝色、绿色的占比。这样的话，由于 color 类需要实现的功能和三维向量类高度相似，就能直接 using Color = Vector3<double> 了。最好，再实现一个把 color 类输出的函数。

假设将三维向量类和颜色类分别封装在 Vector3.hpp 和 Color.hpp 中，那么可以改写第一部分中我们的第一个程序为一个更工程的实现：

#include "Vector3.hpp"
#include "Color.hpp"

int main() {
    int n = 256, m = 256;
    std::cout << "P3\n" << n << " " << m << "\n255\n";
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j) {
            Color c = Color(0, double(i)/(n-1), double(j)/(m-1));
            printColor(std::cout, c);
        }
    return 0;
}

这个程序的运行结果应当要和上面例子里的图片一样。

虚拟视口

通常情况下，我们并不只需要 $256\times 256$ 的图片，而且，也不是只是需要 $1:1$ 大小的图片，所以我们需要设置好图片的宽高比。例如 $16:9$.

实现方法很简单，我们可以设置一个宽高比和图片的宽度（或高度），然后计算图片的高度（或宽度）。需要避免非常极端的数据，导致图片另一维没有数据（如高度为零）。

但由于我们要用（最好用）int 去存图片的大小，所以这时候我们设置的宽高比 $16:9$ 并不一定是真实的宽高比。所以，在需要用到图片宽高比的计算时，请还是使用宽度除以高度去计算，而不是直接使用预设值。

我们希望构建一个假想的“屏幕”，作为一个虚拟视口，在这个屏幕上显示我们能看到的东西：

那么我们需要设置这个屏幕的期望宽高比，以及这个屏幕的大小。在后续的代码里，我设定的宽高比为 $16:9$，画面像素大小为 $800\times 450$.

为了方便应对不同大小的情况，最好不要在程序的各个地方散落屏幕的实际大小 $800\times 450$，不然改不动。那么好的方法是，我们统一一下这个大小，将屏幕的高度统一为 $2$，对应坐标 $[-1,1]$ 的值（即向量单位化之后，任意一维的取值范围），宽度依然用高度和宽高比计算。

接下来我们需要一个摄像机，也就是图片中的演员。显然，为了方便，将摄像机设置为空间的原点。为了符合数学和物理直觉，我们建右手系（Unity 里面建的是左手系，那很抽象了）。为了一些原因，我们设置 $y$ 轴向上、$x$ 轴向右，$z$ 轴的负方向指向虚拟视口。此外，定义摄像机到虚拟视口的长度为焦距，这个值我们设为 $1$.

架好摄像机

为了能让我们看到东西，上帝说，要有光。

一束光线，在正常理解中应该是一条射线。用 $R(t): A+t\vec v$ 表示，其中 $t\geqslant 0$。其中 $A$ 是起始点，$\vec v$ 是方向向量。这个东西封装起来很显然。为了方便，我额外增加了一个函数，用于返回射线在 $t$ 取一个值的时候的结果。

下面我们考虑从一个固定点（摄像机）出发，发出的光线（视线），因为我们已经有了虚拟视口的宽度和高度了，那么对于每个像素点，我们都可以引出一条光线（下图只画了五条，仅作示例）：

对于每条光线，我们其实在考虑下面的问题：找到这条光线和空间中物体的交点，并提取第一个交点处的颜色信息。把这个颜色信息作为对应像素点的颜色呈现出来。

那么到目前为止，可以写出下面的程序。我们期望使用函数 getRayColor 来实现找到一条光线和空间中和物体第一次相交的交点的颜色。不过我们可以先画画简单的东西试试，比如一个渐变（程序中的 _White, _Blue 是我定义的常颜色）：

#include "Color.hpp"
#include "Ray.hpp"
#include "Vector3.hpp"

Color getRayColor(const Ray &ray) {
    // simple example: radial gradient
    double k = ray.dir.x * ray.dir.x + ray.dir.y * ray.dir.y;
    return (1 - k) * _White + k * _Blue;
}

int main() {
    // Viewport
    double expectRatio = 16.0 / 9.0;
    int width = 800;
    int height = std::max(1, (int)(width / expectRatio));
    double viewHeight = 2.0;
    double viewWidth = (double)(width) / height * viewHeight;

    // Camera
    double focalLengh = 1.0;
    Vec3 cameraPos(0, 0, 0);

    Vec3 deltaX = Vec3(viewWidth / width, 0, 0);
    Vec3 deltaY = Vec3(0, -viewHeight / height, 0);
    Vec3 pixUpLeftPos = cameraPos - Vec3(viewWidth / 2, -viewHeight / 2, focalLengh) + 0.5 * (deltaX + deltaY);

    // Render
    std::cout << "P3\n" << width << " " << height << "\n255\n";
    for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++j) {
            Vec3 pixPos = pixUpLeftPos + deltaX * j + deltaY * i;
            Vec3 rayDir = (pixPos - cameraPos).unit();
            Ray ray(cameraPos, rayDir);
            Color col = getRayColor(ray);
            printColor(std::cout, col);
        }

    return 0;
}

大概的逻辑就是先设置好视口，然后架一台摄像机，为了方便，计算相邻两个像素点之间的位移和整个视口的左上角的像素的位移，这样就可以根据当前的视口坐标 $(i,j)$ 确定这个像素的中心点坐标 picPos，从而确定光线。

对于得到的图案，可以先理论分析一下，由于我们单位化了光线的方向向量，所以 $x^2+y^2+z^2=1$. 我们取 $x^2+y^2$ 为一个渐变参数，极限情况为 $z=1$（即光线垂直打到虚拟视口平面上）和 $z=0$（光线和视口平面平行，这不可能发生）。此时对应的颜色为白色和蓝色。因此最终图案是一个从中间到四周，从白色到蓝色的辐射渐变。

实际效果也是如此：

简单的实例：球

下面我们开始真的考虑在坐标系里面加点东西。出于简单，先试试加一个球体。

我们都知道球心在 $C(x_c,y_c,z_c)$，半径为 $r$ 的球的公式是 $(x-x_c)^2+(y-y_c)^2+(z-z_c)^2=r^2$. 但现在，我们希望知道光线 $R: A+t\vec v$ 和这个球有没有交点，如果有，还要计算出靠近 $A$ 这一侧的交点。假设个交点是 $P=A+t\vec v$，则：

$\begin{aligned} & (C-(A+t\vec v))\cdot(C-(A+t\vec v)) = r^2 \\ \Rightarrow& t^2\vec v\cdot\vec v - 2t\vec v\cdot(C-A) + (C-A)\cdot(C-A) - r^2 = 0 \\ \Rightarrow& \begin{cases} t = \frac{-b\pm\sqrt{b^2-4ac}}{2a} \\ a = \vec v\cdot\vec v \\ b = -2\vec v\cdot(C-A) \\ c = (C-A)\cdot(C-A)-r^2 \end{cases} \end{aligned}$

实际上，这个结果可以稍微化简一下，卡个常数：

$\begin{aligned} & \begin{cases} a = \vec v\cdot\vec v \\ h:= \vec v\cdot(C-A) \\ b = -2h \\ c = (C-A)\cdot(C-A)-r^2 \end{cases} \\ \Rightarrow t &= \frac{-b\pm\sqrt{b^2-4ac}}{2a} \\ &= \frac{2h\pm\sqrt{4h^2-4ac}}{2a} \\ &= \frac{h\pm\sqrt{h^2-ac}}{a} \end{aligned}$

于是我们可以针对一个球，计算光线是否和球相交，若相交，还能得出一个合法的 $t$ 的值。

由于后续可能会加入更多的球，或其他什么乱七八糟的形状，我们充分发挥 C++ 的多态，设置一个父类（基类）来统称这些玩意。最开始我是想把类型命名为 Object 的，但是考虑到实际游戏里有很多东西不能被光线探测到，比如空气和碰撞箱，所以不会命名了。然而《Ray Tracing in One Weekend》这本书里面给出了一个很妙的方式：为了体现这些物体是会被光线击中的这一特性，命名为 Hittable.

在父类 Hittable 下定义子类球体，即 Sphere：

class Hittable {
  public:
    virtual double hit(const Ray &ray) const = 0;
};

class Sphere : public Hittable {
  public:
    Vec3 C;
    double r;
    Sphere(const Vec3 &_C = Vec3(0, 0, 0), double _r = 0) : C(_C), r(_r) {}

    double hit(const Ray &ray) const override {
        double a = dot(ray.dir, ray.dir);
        double h = dot(ray.dir, C - ray.A);
        double c = dot(C - ray.A, C - ray.A) - r * r;
        double delta = h * h - a * c;
        if(delta < 0)
            return -1;
        double sqDelta = std::sqrt(delta);
        if (h + sqDelta < 0)
            return -1;
        return (h <= sqDelta ? (h + sqDelta) / a : (h - sqDelta) / a);
    }
};

现在的球没有任何颜色，非常的单调，我们可以想一个好的方法给球上一个花里胡哨的颜色。例如，我们取光线和球的交点处在球表面上的向外的法向量并单位化，那么单位化后的法向量，$x,y,z$ 分别属于 $[-1,1]$，可以设计映射去和 $r,g,b$ 去一一对应。这样子，一个点的颜色和法向量一一对应了，虽然说不出哪里有意思但就是很有意思就是了。

关于法向量的计算，显然，法向量的方向就是 $P-C$ 的方向.

这里引入法向量，因为想想也觉得是，我们后面肯定会需要用到交点在面上的法向量去算一些东西的。

底色（即没有击中物体的光线的颜色）直接采用你喜欢的渐变色即可（请注意，这个程序实现的渐变不是垂直渐变）：

Sphere S(Vec3(0, 0, -2), 0.5);

Color getRayColor(const Ray &ray) {
    double t = S.hit(ray);
    if (t >= 0){
        Vec3 n = (ray.at(t) - S.C).unit();
        return 0.5 * Color(n.x + 1, n.y + 1, n.z + 1);
    }
    double k = 0.5 * (ray.dir.y + 1);
    return (1 - k) * _White + k * Color(0.7, 0.7, 1);
}

那么得到的图案应该是：

实例集合封装

为了方便我们实现一大堆的实例的碰撞检测，我们可以考虑封装一个实例集合。为了方便不区分地直接使用 hit 去判定一个光线是否和这个这个集合内地某个实例相交，我们干一个可能反人类地操作：把这个集合视为 Hittable 的一个子类。

using std::make_shared;
using std::shared_ptr;

class HittaleSet : public Hittable {
  public:
    std::vector<shared_ptr<Hittable>> vec;
    HittaleSet() {}

    void clear() {
        vec.clear();
    }
    void insert(shared_ptr<Hittable> obj) {
        vec.push_back(obj);
    }
    double hit(const Ray &ray) const override {
        double t = -1;
        for (const auto &x : vec) {
            double res = x->hit(ray);
            if (res >= 0 && (t < 0 || res < t))
                t = res;
        }
        return t;
    }
};

注意

由于 Hittable 是一个父类（基类），如果你往里面加入子类对象的时候，会触发切片，导致只有父类自己的内容会被保留，子类自定义的东西就没有了。你的子类的 hit 函数就会原地蒸发。

因此，我们用 shared_ptr，不能直接使用 vector<Hittable>.

摄像头封装

同样的，考虑到很可能不仅仅使用一个固定在原点的摄像头，我们需要封装一个摄像头类。这一段很轻松，其实就是把主函数的东西移植到类里面。

……想法很美好，但我们很快发现主函数里面的有一个东西是难以改造的：

Sphere S(Vec3(0, 0, -2), 0.5);

Color getRayColor(const Ray &ray) {
    double t = S.hit(ray);
    if (t >= 0){
        Vec3 n = (ray.at(t) - S.C).unit();
        return 0.5 * Color(n.x + 1, n.y + 1, n.z + 1);
    }
    double k = 0.5 * (ray.dir.y + 1);
    return (1 - k) * _White + k * Color(0.7, 0.7, 1);
}

这个部分，尤其是计算法向量的部分，是针对了对应的物体是一个球计算的。但实际上，我们需要支持更多类型的实例时，就不能单纯这么算了。球有球的法向量算法，正方体有正方体的，各有各的，而且传入的类型的 Hittable，所以怎么算法向量还得额外判断对应的子类。

这实在是太麻烦了，所以我们需要附带改进一下 hit 函数，使得其同时返回切点和对应的法向量。一种经典而且常数更小的方法是，自定义一个返回值类型，比如 hitResult，给 hit 函数传入一个 hitResult 类型的变量的地址，去把结果存到这个变量里。当然也可以直接返回 hitResult 类型，但是常数更大。因此这里采用前者。

hitResult 类型里，目前只需要两个变量，一个 $t$ 一个法向量 $\vec n$ 就可以了，为了方便，默认值设为 $t=-1$. 然后我们修改 hit 函数：

// Hittable
bool hit(const Ray &ray, hitResult &res) const override {
    hitResult tmp;
    bool ret = false;
    for (const auto &x : vec)
        if(x->hit(ray, tmp) && (res.t < 0 || tmp.t < res.t)){
            res = tmp;
            ret = true;
        }
    return ret;
}
// Sphere
bool hit(const Ray &ray, hitResult &res) const override {
    double a = dot(ray.dir, ray.dir);
    double h = dot(ray.dir, C - ray.A);
    double c = dot(C - ray.A, C - ray.A) - r * r;
    double delta = h * h - a * c;
    if (delta < 0)
        return false;
    double sqDelta = std::sqrt(delta);
    if (h + sqDelta < 0)
        return false;
    res.t = h <= sqDelta ? (h + sqDelta) / a : (h - sqDelta) / a;
    res.n = (ray.at(res.t) - C).unit();
    return true;
}

这样我们就可以实现摄像机类了：

class Camera {
  private:
    double expectRatio;
    int width, height;
    double viewHeight, viewWidth;
    double focalLength;
    Vec3 pos;
    Vec3 deltaX, deltaY, pixUpLeftPos;

    Color getRayColor(const Ray &ray, const Hittable &objs) {
        hitResult res;
        if (objs.hit(ray, res))
            return 0.5 * Color(res.n.x + 1, res.n.y + 1, res.n.z + 1);
        double k = 0.5 * (ray.dir.y + 1);
        return (1 - k) * _White + k * Color(0.7, 0.7, 1);
    }

  public:
    Camera(double _expectRatio = 16.0 / 9.0, int _width = 800, double _viewHeight = 2.0, double _focalLength = 1.0,
           Vec3 _pos = Vec3(0, 0, 0))
        : expectRatio(_expectRatio), width(_width), height(std::max(1, (int)(_width / _expectRatio))),
          viewHeight(_viewHeight), viewWidth((double)(_width) / height * viewHeight), focalLength(_focalLength),
          pos(_pos), deltaX(viewWidth / width, 0, 0), deltaY(0, -viewHeight / height, 0),
          pixUpLeftPos(_pos - Vec3(viewWidth / 2, -viewHeight / 2, focalLength) + 0.5 * (deltaX + deltaY)) {}

    void render(const Hittable &objs) {
        std::cout << "P3\n" << width << " " << height << "\n255\n";
        for (int i = 0; i < height; ++i)
            for (int j = 0; j < width; ++j) {
                Vec3 pixPos = pixUpLeftPos + deltaX * j + deltaY * i;
                Vec3 rayDir = (pixPos - pos).unit();
                Ray ray(pos, rayDir);
                Color col = getRayColor(ray, objs);
                printColor(std::cout, col);
            }
    }
};

随后我们的主函数大大简化：

int main() {
    HittaleSet objs;
    objs.insert(make_shared<Sphere>(Vec3(0, 0, -2), 0.5));
    Camera cam(16.0 / 9.0, 800, 2.0, 1.0, Vec3(0, 0, 0));
    cam.render(objs);
    return 0;
}

运行主函数，应该会得到和前面一样的彩色球的图案。

抗锯齿

即使你不放大那个彩球图象也会发现，那张图的锯齿化非常严重。这是因为每个像素要么对应的光线打到球上了，此时整个像素的颜色就是那一条打到球上的光线的颜色，要么没有打到，取背景色。用一条光线独断导致边界像素突出。

那么就有了一个抗锯齿的大致的方向：就是让像素的颜色稍微“中和”一下，使得边界处的颜色过渡稍微平缓。

注意到对于每个像素方格，我们只取了格子的正中心的一个点，读取这个点对应的光线的颜色。结合上面的思想，我们可以换成采样：在这个格子（或比这个格子稍大一点的范围，设置合理即可）内额外采样一些点去得到一个，把这些颜色平均起来。这样子一个像素点对应的颜色期望下应该是这个像素对应的实际范围的一个“平均”颜色。

为此，我们需要给摄像机加入一个新的参数 nSample 表示单个像素的额外采样量。

为了实现随机，我新建了一个头文件 Math.hpp，将和数学有关的操作和常量（包括之前的计算几何全家桶）放进去。在里面我实现了随机数生成器和辅助函数，由于实现简单，此处忽略。

下面是修改后的渲染，我设置的范围比这个像素格子的大小稍大一点：

Color getSampleColor(const Vec3 &pixPos, const Hittable &objs) {
    Vec3 samPos = pixPos + deltaX * rnd.rand<double>(-0.7, 0.7) + deltaY * rnd.rand<double>(-0.7, 0.7);
    Ray ray(pos, (samPos - pos).unit());
    return getRayColor(ray, objs);
}

void render(const Hittable &objs) {
    double sampleRatio = 1.0 / (nSample + 1);
    std::cout << "P3\n" << width << " " << height << "\n255\n";
    for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++j) {
            Vec3 pixPos = pixUpLeftPos + deltaX * j + deltaY * i;
            Ray ray(pos, (pixPos - pos).unit());
            Color col = getRayColor(ray, objs);
            for(int k = 0; k < nSample; ++k)
                col += getSampleColor(pixPos, objs);
            col *= sampleRatio;
            printColor(std::cout, col);
        }
}

将采样数设置为 $50$ 时，运行修改后的渲染，可以得到这样的图片：

自行放大和上一张彩色球图片对比，可以发现抗锯齿发挥作用了。

光追

材料类

我们希望能实现不同材质的球，有些只是单纯的球（触发漫反射），有些可能是镜面球，有些可能是玻璃球……球的种类影响光线的逻辑。所以我们需要一个材料类来管理不同的材料。

一个材料的信息应该包含：如果给定入射光线、法向量等信息（hitResult），会不会产生反射光线，如果会，会产生怎样的反射光线，又吸收了多少光线。

class Material {
  public:
    virtual bool scatter(const Ray &rayIn, const hitResult &res, Color &attn, Ray &rayOut) const = 0;
};

hitResult 中需要添加光线碰到的实例的材料信息，但是材料类里面也用到了 hitResult，为了避免头文件互相引用或者其他冲突，在 hitResult 对应的头文件中，只声明材料类，不引用材料类头文件即可。

然后我们修改 getRayColor 的逻辑，为了避免大量的反射导致超绝递归，我们还要给摄像头类设置一个最大递归深度，如果超过这个递归深度，则视为没有光（黑色），调用 getRayColor 时，dep 的值设置为最大递归深度：

Color getRayColor(const Ray &ray, const Hittable &objs, int dep) {
    if (dep == 0)
        return _Black;
    hitResult res;
    if (objs.hit(ray, res)) {
        Ray scat;
        Color attn;
        if (res.mat->scatter(ray, res, attn, scat))
            return atten * getRayColor(scat, objs, dep - 1);
        return _Black;
    }
    double k = 0.5 * (ray.dir.y + 1);
    return (1 - k) * _White + k * Color(0.7, 0.7, 1);
}

别忘记在球的类型中，添加球的材质。

漫反射（郎伯模型）

现在我们实现一个会进行漫反射的材料。当光线照射到漫反射材料的时候，我们用郎伯模型去计算反射情况：

假设物体表面平滑，入射光完全漫反射后在所有方向的辐射亮度相同，漫反射的光强和入射角的余弦成正比，和反射方向无关。

需要注意的是，虽然我们代码实现是从摄像头发出光线的，但事实显然并非如此，是光线反射到摄像头的。请不要弄反入射角和反射角。

我们可以用如下方法轻松地构造和入射角（从摄像头的视角，或许是“反射角”？）的余弦成正比分布的随机光线：

1. 得到交点 $P$ 和法向量（单位化过的）$\vec n$ 后，得到点 $O=P+\vec n$. 以 $O$ 为球心作单位球。
2. 取单位球上随机一个不为 $P$ 的点 $S$.
3. 连接 $PS$，$PS$ 就是我们需要的光线。

以上是示意图。

这种随机方法符合郎伯分布的证明

以 $\vec n$ 所在轴为极轴建立极坐标系，假设 $S$ 的和极轴夹角为 $\varphi$（在上图中已经标记），方位角为 $\phi$. 那么固定 $\varphi$ 得到的一条“纬线”的环面面积

$\mathrm{d} A = 2\pi\sin\varphi\,\mathrm{d}\varphi$

球的表面积为 $4\pi$，于是 $\varphi_0$ 落在 $[\varphi,\varphi+\mathrm{d}\varphi]$ 的概率为

$\frac{\mathrm{d}A}{4\pi} = \frac{\sin\varphi}{2}\,\mathrm{d}\varphi = \sin(2\theta)\,\mathrm{d}\theta$

对于给定 $\theta$ 和任意方位角 $\phi$，有固体角面积

$\mathrm{d}\Omega = 2\pi\sin\theta\,\mathrm{d}\theta$

于是得到概率密度

$p = \frac{\sin(2\theta)\,\mathrm{d}\theta}{\mathrm{d}\Omega} = \frac{\cos\theta}{\pi}\propto \cos\theta$

显然，漫反射本身还要吸收一些光，在这里我们先随便设置一个默认吸收量，例如 $50\%$. 还有一个需要注意的点：我们需要判断当前的光线到底照射到的是球的内表面还是外表面，所以需要判断法向量和光线的同向关系。又还有一个需要注意的点，我们不希望随机取到的点 $S$ 离 $P$ 太近，造成一系列 bug，所以需要特判。

根据以上理论，建立发生漫反射的材料 Diffuse 类：

class Diffuse : public Material {
  private:
    Color attenCol;

  public:
    Diffuse() : attenCol(1, 1, 1) {}
    Diffuse(const Color &_attenCol) : attenCol(_attenCol) {}
    bool scatter(const Ray &rayIn, const hitResult &res, Color &atten, Ray &rayOut) const override {
        Vec3 realN = dot(rayIn.dir, res.n) > 0 ? -res.n : res.n;
        Vec3 outDir = (realN + rnd.randUnit()).unit();
        if (outDir.nearZero())
            outDir = realN;
        rayOut = Ray(rayIn.at(res.t), outDir);
        atten = attenCol;
        return true;
    }
};

由于我们目前还只有一个球，为了方便检查，我们增加一个巨大的球，用来假装是地面。我设置小球为 $C:(0,0,-2),r=0.5$，地面球为 $C:(0,-100.5,-2),r=100$. 两个球的光线反射率为 $70\%$，即吸收量为 $30\%$.

可以勉强看到接触部分有更深的阴影。

阴影失真

为什么这个东西明明只有 $30\%$ 的吸收率，但是还是这么黑？

《Ray Tracing in One Weekend》提到了这样的阴影失真的问题。这是一个非常小但后果很严重的错误：我们使用 at 试图得到光线和球体表面的交点 $P$ 的精确坐标，但是由于浮点数误差，我们算出来的 $P$ 可能发生了些微偏移。这个偏移有可能刚好从球的表面偏移到了内部，这样的话新的光线会从很贴近边界的球的一个内点发射，然后马上和球的表面形成二次相交。

由于浮点数误差难以避免，所以一个比较好的解决这个问题的方案是：直接忽略和出发点非常接近的交点。

我们定义区间类 Interval，包含 $l,r$ 两个参数，表示一个闭区间 $[l,r]$，用于方便修改程序。然后修改 hit 函数和 getRayColor 函数：

// Hittable
bool hit(const Ray &ray, Interval rangeT, hitResult &res) const override {
    hitResult tmp;
    double nowR = rangeT.r;
    for (const auto &x : vec)
        if (x->hit(ray, Interval(rangeT.l, nowR), tmp) && (res.t < 0 || tmp.t < res.t)) {
            nowR = tmp.t;
            res = tmp;
        }
    return rangeT.contain(res.t);
}
// Sphere
bool hit(const Ray &ray, Interval rangeT, hitResult &res) const override {
    double a = dot(ray.dir, ray.dir);
    double h = dot(ray.dir, C - ray.A);
    double c = dot(C - ray.A, C - ray.A) - r * r;
    double delta = h * h - a * c;
    if (delta < 0)
        return false;
    double sqDelta = std::sqrt(delta);
    double t1 = (h - sqDelta) / a, t2 = (h + sqDelta) / a;
    if (rangeT.contain(t1))
        res.t = t1;
    else if (rangeT.contain(t2))
        res.t = t2;
    else
        return false;
    res.n = (ray.at(res.t) - C).unit();
    res.mat = mat;
    return true;
}

getRayColor 函数只是在 hit 判断中多传入一个区间参数 $[\varepsilon_0,+\infty]$，我的程序中取了 $\varepsilon_0=10^{-4}$. 渲染结果如下：

可以发现，球过黑问题有明显改善。

关于为什么我们的球没有设置颜色，但是它渲染出来看起来是淡紫色的，因为我的例子中，背景板是紫色的。

Gamma 校正

虽然这个球已经很亮了，但那是对比的阴影失真的版本。实际上，反射率 $70\%$ 的球应该比这更亮。这是因为我们目前的图像没有进行 Gamma 校正。

存在一些客观因素导致人眼看到的亮度和计算机实际意图显示的亮度并不是线性对等关系。即你给它输入 $x\,(0\leqslant x\leqslant 1)$，它输出一个亮度 $y\,(0\leqslant y\leqslant 1)$，其中 $y=x$（即线性空间），则在人眼看来，这个亮度其实不到 $y$（变暗了）。我们感觉的亮度更符合一个非线性函数 $y_1=x^\gamma$，其中 $\gamma$ 通常取 $2.2$（一些科学研究的结果）。为了校正这一点，我们输入亮度 $x$ 的时候，计算机视为意图输出的亮度是 $y_2=x^{\frac{1}{\gamma}}$（即 Gamma 空间），这样就补正了结果。

我们目前输出图像是线性空间内的，所以渲染出来的东西其实在数据上是对的，但是肉眼看起来上和数据相比偏暗。因此我们需要进行 Gamma 补正。这非常好实现，在此忽略修改代码。

我们会得到反射率为 $70\%$ 时，这种东西真正的样子：

镜面反射

镜面反射的逻辑比漫反射好实现得多。那为什么先说漫反射呢，废话，一个全是镜子的图我还看个啥，看不动啊。

由显然的公式，当然画个图也可以知道，已知入射光线方向向量 $\vec v$ 和法向量 $\vec n$（都已经单位化），则反射光线的方向向量为 $\vec v - 2(\vec v\cdot\vec n)\vec n$.

为了方便起见，我们定义的镜子材料还是可能存在颜色吸收的。默认为全反射。

class Mirror : public Material {
  private:
    Color attenCol;

  public:
    Mirror() : attenCol(1, 1, 1) {}
    Mirror(const Color &_attenCol) : attenCol(_attenCol) {}
    bool scatter(const Ray &rayIn, const hitResult &res, Color &atten, Ray &rayOut) const override {
        Vec3 realN = dot(rayIn.dir, res.n) > 0 ? -res.n : res.n;
        Vec3 outDir = rayIn.dir - 2 * dot(rayIn.dir, realN) * realN;
        rayOut = Ray(rayIn.at(res.t), outDir);
        atten = attenCol;
        return true;
    }
};

我们使用复杂一点的球来测试程序，例如：

int main() {
    HittaleSet objs;
    objs.insert(make_shared<Sphere>(make_shared<Diffuse>(Color(0.5, 0.6, 0.7)), Vec3(0, 0.1, -2), 0.5));
    objs.insert(make_shared<Sphere>(make_shared<Diffuse>(Color(0.3, 0.3, 0.3)), Vec3(0, -100.5, -2), 100));
    objs.insert(make_shared<Sphere>(make_shared<Mirror>(), Vec3(-1, 0.5, -2), 0.4));
    objs.insert(make_shared<Sphere>(make_shared<Mirror>(Color(0.9, 0.5, 0.3)), Vec3(1, 0, -1.5), 0.3));
    Camera cam(16.0 / 9.0, 800, 2.0, 1.0, Vec3(0, 0, 0), 50, 10);
    cam.render(objs);
    return 0;
}

第一个球是一个有颜色的漫反射球，第二个球为地面，第三个球为全反射镜子球，第四个球为有光线吸收的带颜色的镜子球。渲染的结果如下：

铜镜：模糊反射

金属的表面通常并非和现代的镜子一样可以进行完全的镜面反射，而是模糊的。

为了实现模糊效果，我们可以不把反射光线的方向唯一确定地算出来，而是在原来的基础上加上一个随机扰动，这样反射光线就会在一个区域内，从而实现平均效果，也就实现了模糊效果。

由于模糊反射的代码实现和镜面反射高度相似，我们可以考虑把它们结合起来。我们把 Mirror 类重命名为反光物体类 Glossy，加入一个新的参数 $fuzz$，表示随机扰动的幅度。特别地，$fuzz=0$ 为不扰动，此时和正常的镜面反射无异。

class Glossy : public Material {
  private:
    Color attenCol;
    double fuzz;

  public:
    Glossy(double _fuzz = 0) : attenCol(1, 1, 1), fuzz(_fuzz) {}
    Glossy(const Color &_attenCol, double _fuzz = 0) : attenCol(_attenCol), fuzz(_fuzz) {}
    bool scatter(const Ray &rayIn, const hitResult &res, Color &atten, Ray &rayOut) const override {
        Vec3 realN = dot(rayIn.dir, res.n) > 0 ? -res.n : res.n;
        Vec3 outDir = rayIn.dir - 2 * dot(rayIn.dir, realN) * realN + fuzz * rnd.randUnit();
        if (dot(outDir, realN) < 0)
            return false;
        rayOut = Ray(rayIn.at(res.t), outDir.unit());
        atten = attenCol;
        return true;
    }
};

我们试试把上述例子的有色镜子球改成有色的 $fuzz=0.3$ 的铜镜球，得到结果：

玻璃：折射

清澈透明的材料会发生光的折射，例如玻璃。当然其实是可能会同时发生折射和反射，也可能只有反射（全反射）。以上都是初中物理内容，只在必要时作提及。

这样的材料有一个折射率。光从折射率 $\eta$ 的介质入射到折射率 $\eta'$ 的介质时，设入射角为 $\theta$，折射角为 $\theta'$，则：

$\frac{\eta}{\eta'} = \frac{\sin\theta'}{\sin\theta}$

我们希望用向量的方式算出来，现在有入射光线的方向向量 $\vec v$ 和法向量 $\vec n$，则折射光线的方向向量 $\vec d$：

$\begin{aligned} & \eta(\vec v\times \vec n) = \eta'(\vec d\times \vec n) \\ \Rightarrow & (\eta\vec v - \eta'\vec d)\times\vec n = \bf 0 \end{aligned}$

这说明 $\eta\vec v - \eta'\vec d$ 和法向量共线. 说明其在垂直法向量方向的分向量 $\eta\vec v_x = \eta'\vec d_x$，而 $\vec v_x = \vec v-(\vec v\cdot\vec n)\vec n$. 于是：

$\begin{aligned} \vec d_x &= \frac{\eta}{\eta'}[\vec v-(\vec v\cdot\vec n)\vec n] \\ \vec d_y &= -\vec n\sqrt{1-|\vec d_x|^2} \\ \vec d &= \vec d_x + \vec d_y \end{aligned}$

于是我们可以计算出折射光线。

但实际上，玻璃还会反射。正常情况下，反射率 $R$ 会随着角度变化。这个叫菲涅尔效应，这是一个非常恐怖的方程（可以阅读这篇知乎：菲涅尔方程（Fresnel Equation） - 知乎）。好在我们还有一个神仙，叫 Schlick 近似，这使得我们可以用一个多项式函数来近似计算菲涅尔方程。

首先我们需要一个基础反射率 $R_0$，这是入射角为 $0$ 时的反射率，根据菲涅尔方程有：

$R_0 = \left( \frac{\eta-\eta'}{\eta+\eta'} \right)^2$

Schlick 近似的公式是，对于入射角 $\theta$，有：

$R_\theta = R_0 + (1-R_0)(1-\cos\theta)^5$

类似地，我们可能需要毛玻璃的效果，即磨砂玻璃，所以和前面的模糊反射一样，我们设置 $fuzz$ 来进行模糊处理。可以考虑只写折射光模糊，也可以同时写折射模糊和反射模糊，我这里选择两个都写。

除了实心的玻璃球之外，我们还可以考虑手动实现空心的玻璃球，但这是没必要的。单纯使用两个球体就可以解决：假如定义一个折射率为 $\eta_0$ 的玻璃球，那么取一个合适的厚度，在玻璃球里面定义一个假装是空气球的玻璃球即可。空气球的折射率应该为 $\frac{1}{\eta_0}$.

下面是玻璃类的实现：

class Glass : public Material {
  private:
    Color attenCol;
    double eta;
    double R0;
    double fuzzReflect, fuzzRefract;
    double getR(const double &cosTheta) const {
        return R0 + (1 - R0) * pow(1 - cosTheta, 5);
    }

  public:
    Glass(double _eta = 1.5, double _fuzzReflect = 0, double _fuzzRefract = 0)
        : attenCol(1, 1, 1), eta(_eta), R0((_eta - 1) * (_eta - 1) / (_eta + 1) / (_eta + 1)),
          fuzzReflect(_fuzzReflect), fuzzRefract(_fuzzRefract) {}
    Glass(const Color &_attenCol, double _eta = 1.5, double _fuzzReflect = 0, double _fuzzRefract = 0)
        : attenCol(_attenCol), eta(_eta), R0((_eta - 1) * (_eta - 1) / (_eta + 1) / (_eta + 1)),
          fuzzReflect(_fuzzReflect), fuzzRefract(_fuzzRefract) {}
    bool scatter(const Ray &rayIn, const hitResult &res, Color &atten, Ray &rayOut) const override {
        double realETA = eta;
        Vec3 realN = res.n;
        if (dot(rayIn.dir, res.n) > 0) {
            realETA = 1.0 / realETA;
            realN = -realN;
        }
        double cosTheta = dot(-rayIn.dir, realN);
        double sinTheta = std::sqrt(1 - cosTheta * cosTheta);
        if (sinTheta * realETA > 1.0 || rnd.rand<double>(0, 1) < getR(cosTheta)) { // reflect
            Vec3 outDir = rayIn.dir - 2 * dot(rayIn.dir, realN) * realN + fuzzReflect * rnd.randUnit();
            if (dot(outDir, realN) < 0)
                return false;
            rayOut = Ray(rayIn.at(res.t), outDir.unit());
            atten = attenCol;
            return true;
        }
        // refract
        Vec3 outX = realETA * (rayIn.dir + realN * cosTheta);
        Vec3 outY = -realN * std::sqrt(std::fabs(1 - outX.len2()));
        Vec3 outDir = outX + outY + fuzzRefract * rnd.randUnit();
        if (dot(outDir, realN) > 0)
            return false;
        rayOut = Ray(rayIn.at(res.t), outDir.unit());
        atten = attenCol;
        return true;
    }
};

我们在主函数中添加两个玻璃类型的球：

int main() {
    HittaleSet objs;
    objs.insert(make_shared<Sphere>(make_shared<Diffuse>(Color(0.5, 0.6, 0.7)), Vec3(0, 0.1, -2), 0.5));
    objs.insert(make_shared<Sphere>(make_shared<Diffuse>(Color(0.3, 0.3, 0.3)), Vec3(0, -100.5, -2), 100));
    objs.insert(make_shared<Sphere>(make_shared<Glossy>(), Vec3(-1, 0.5, -2), 0.4));
    objs.insert(make_shared<Sphere>(make_shared<Glossy>(Color(0.9, 0.5, 0.3), 0.3), Vec3(1, 0, -1.5), 0.3));
    objs.insert(make_shared<Sphere>(make_shared<Glass>(), Vec3(0.7, -0.2, -1), 0.3));
    objs.insert(make_shared<Sphere>(make_shared<Glass>(Color(0.7, 0.9, 0.8), 1.5, 0.3, 0.3), Vec3(-0.5, -0.1, -1), 0.2));
    Camera cam(16.0 / 9.0, 800, 2.0, 1.0, Vec3(0, 0, 0), 50, 10);
    cam.render(objs);
    return 0;
}

渲染结果为：

万能的区域类和程序优化

上面的程序跑出来的这个玻璃球看起来挺怪的。尤其是左边的有色玻璃，不管是长相还是颜色。对于长相，像一个鸡蛋，里面一圈深色，外面颜色较浅。对于颜色，实际上，这个颜色是错误的。我们先来考虑简单一点的，即颜色错误的问题。这是因为我们现在的颜色衰减，是根据反射或折射的次数决定的，但现实里，显然是由光线在介质中传播的距离决定的。

为了实现更真实的有色玻璃，我们需要想个办法能计算光线在这个几何体里面走了多远。——而且最好能避免直接把材料类和几何类直接耦合起来，那写死人了。

为了实现这一点，我们写一个额外的区域类。这样就能方便我们实现光线在介质中传播的距离，同时也方便实现后续实现一些“飘渺”的东西，比如体积雾。

区域类负责干的事情是处理光线在体积内的移动距离。当然我们顺便把颜色吸收的系数也可以移植进来，更加方便计算。

对于一个区域，我们希望能记录这个区域的边界（几何类）、材质（材料类）和颜色吸收量。我们希望用光线的路径长度 $dis$ 和单位衰减量 $att_0$ 来计算总的颜色衰减 $atten$，则根据 Beer-Lambert 定律，有：

$atten = e^{-dis\cdot att_0}$

需要注意的是，这是光线能穿过区域的情况，即区域透明，显然不能用在漫反射和镜面反射材料上。这是很好解决的问题：给每个材料，标记一下是否透明即可。

class wrapZone : public Hittable {
  public:
    shared_ptr<Hittable> bnd;
    shared_ptr<Material> mat;
    Color att0;
    wrapZone(shared_ptr<Hittable> _bnd, shared_ptr<Material> _mat, const Color &_col)
        : bnd(_bnd), mat(_mat), att0(_mat->isTrans ? _White - _col : _col) {}
    bool hit(const Ray &ray, Interval rangeT, hitResult &res) const override {
        return bnd->hit(ray, rangeT, res);
    }
    Color atten(const Ray &ray) const {
        if (!mat->isTrans)
            return att0;
        hitResult res1, res2;
        if (!bnd->hit(ray, Interval(-INFINITY, INFINITY), res1))
            return Color(1, 1, 1);
        if (!bnd->hit(ray, Interval(res1.t + 0.0001, INFINITY), res2))
            return Color(1, 1, 1);
        double dis = res2.t - res1.t;
        return Color(std::exp(-att0.x * dis), std::exp(-att0.y * dis), std::exp(-att0.z * dis));
    }
    bool scatt(const Ray &rayIn, const hitResult &res, Ray &rayOut) const {
        return mat->scatter(rayIn, res, rayOut);
    }
};

由于我们在这里定义了颜色衰减率 $att_0$，材料类的颜色可以进行删去。我们删除材料类的颜色吸收参数之后，将主函数内所有物体的定义改为使用区域类定义，这样的话，每个几何体的材料也可以从几何体类的参数列表里删去。最终主函数得到：

int main() {
    HittaleSet objs;
    objs.insert(make_shared<wrapZone>(
        make_shared<Sphere>(Vec3(0, 0.1, -2), 0.5),
        make_shared<Diffuse>(),
        Color(0.5, 0.6, 0.7)
    ));
    objs.insert(make_shared<wrapZone>(
        make_shared<Sphere>(Vec3(0, -100.5, -2), 100),
        make_shared<Diffuse>(),
        Color(0.3, 0.3, 0.3)
    ));
    objs.insert(make_shared<wrapZone>(
        make_shared<Sphere>(Vec3(-1, 0.5, -2), 0.4),
        make_shared<Glossy>(),
        Color(1, 1, 1)
    ));
    objs.insert(make_shared<wrapZone>(
        make_shared<Sphere>(Vec3(1, 0, -1.5), 0.3),
        make_shared<Glossy>(0.3),
        Color(0.9, 0.5, 0.3)
    ));
    objs.insert(make_shared<wrapZone>(
        make_shared<Sphere>(Vec3(0.7, -0.2, -1), 0.3),
        make_shared<Glass>(),
        Color(1, 1, 1)
    ));
    objs.insert(make_shared<wrapZone>(
        make_shared<Sphere>(Vec3(-0.5, -0.1, -1), 0.2),
        make_shared<Glass>(1.5, 0.3, 0.3),
        Color(0.7, 0.9, 0.8)
    ));
    Camera cam(16.0 / 9.0, 800, 2.0, 1.0, Vec3(0, 0, 0), 50, 10);
    cam.render(objs);
    return 0;
}

HittableSet 中的元素类型，也从父类 Hittable 改成 wrapZone.

最后重构 getRayColor 的逻辑，因为现在颜色吸收不再从材料里读取，而是从区域类里读取。hitResult 中，材料类改为体积类。然后：

Color getRayColor(const Ray &ray, const HittaleSet &objs, int dep) {
    if (dep == 0)
        return _Black;
    hitResult res;
    if (objs.hit(ray, Interval(0.0001, INFINITY), res)){
        Ray scat;
        auto obj = res.zone;
        if(obj->scatt(ray, res, scat))
            return obj->atten(ray) * getRayColor(scat, objs, dep - 1);
        return _Black;
    }
    double k = 0.5 * (ray.dir.y + 1);
    return (1 - k) * _White + k * Color(0.7, 0.7, 1);
}

下面来考虑鸡蛋的形状。我们修改毛玻璃的逻辑，不要在调整后不符合条件就返回，而是找到一个符合条件的结果返回：

class Glass : public Material {
  private:
    double eta;
    double R0;
    double fuzzReflect, fuzzRefract;
    double getR(const double &cosTheta) const {
        return R0 + (1 - R0) * pow(1 - cosTheta, 5);
    }

  public:
    Glass(double _eta = 1.5, double _fuzzReflect = 0, double _fuzzRefract = 0)
        : eta(_eta), R0((_eta - 1) * (_eta - 1) / (_eta + 1) / (_eta + 1)), fuzzReflect(_fuzzReflect),
          fuzzRefract(_fuzzRefract), Material(true) {}
    bool scatter(const Ray &rayIn, const hitResult &res, Ray &rayOut) const override {
        double realETA = eta;
        Vec3 realN = res.n;
        if (dot(rayIn.dir, res.n) > 0) {
            realETA = 1.0 / realETA;
            realN = -realN;
        }
        double cosTheta = dot(-rayIn.dir, realN);
        double sinTheta = std::sqrt(1 - cosTheta * cosTheta);
        if (sinTheta * realETA > 1.0 || rnd.rand<double>(0, 1) < getR(cosTheta)) { // reflect
            Vec3 outDir = rayIn.dir - 2 * dot(rayIn.dir, realN) * realN + fuzzReflect * rnd.randUnit();
            while (dot(outDir, realN) < 0)
                outDir = rayIn.dir - 2 * dot(rayIn.dir, realN) * realN + fuzzReflect * rnd.randUnit();
            rayOut = Ray(rayIn.at(res.t), outDir.unit());
            return true;
        }
        // refract
        Vec3 outX = realETA * (rayIn.dir + realN * cosTheta);
        Vec3 outY = -realN * std::sqrt(std::fabs(1 - outX.len2()));
        Vec3 outDir = outX + outY + fuzzRefract * rnd.randUnit();
        while (dot(outDir, realN) > 0)
            outDir = outX + outY + fuzzRefract * rnd.randUnit();
        rayOut = Ray(rayIn.at(res.t) - 0.001 * realN, outDir.unit());
        return true;
    }
};

运行得到图片：

质量有明显提高。