4 - 下周的光追Part1

发表于 2024-08-25 更新于 2024-09-01 分类于图形学，光追三部曲阅读次数： Waline：本文字数： 14k 阅读时长 ≈ 11 分钟

开始学习大名鼎鼎的光追三部曲系列中的：Ray Tracing: The Next Week！希望我能坚持下去吧。

运动模糊

接下来看看运动模糊。对于真正的相机，快门会在短时间内保持开启，在这段时间内如有物体移动，就会产生运动模糊。

时空(spacetime)光追简介

可以通过在快门打开时的某一瞬射出一条随机光线，这样就能获得一个简单光子的随机估计。只要我们能够决定那一瞬物体的位置，那么就能精准测量出那一瞬的光。这也是随机光线追踪（如蒙特卡洛）最终这么简单的原因。

因为光线追踪器可以确保每一条光线对应的物体位置，那么求交的变化就不大。接下来需要为每条光线存储确切的时间，修改Ray类如下：

class Ray 
{
public:
	Ray() {}
	Ray(const Point3& origin, const Vec3& direction, double time)
		: orig(origin), dir(direction), tm(time) {}
	Ray(const Point3& origin, const Vec3& direction)
		: Ray(origin, direction, 0) {}

	const Point3& origin() const { return orig; }
	const Vec3& direction() const { return dir; }
	double time() const { return tm; }

	// P(t) = Orig + t*dir
	Point3 at(double t) const {
		return orig + t * dir;
	}

private:
	Point3 orig;
	Vec3 dir;
	double tm;
};

管理时间

接下来看看如何管理时间，可以从以下两个角度思考快门的时间：

当前快门的开启时间至下一个快门的开启时间
快门在每帧保持开启的时间

对于标准的电影，通常是24fps；而现代数字电影则由导演决定。每帧都有自己的快门速度，它不需要是每帧的最大持续时间，而是每帧 1/1000 秒，或 1/60 秒。

如果我们想要渲染顺序图，需要给相机设置合适的快门时间：帧间周期，快门持续时间，总帧数。如果相机在移动，场景是静态的那还好；如果场景中有物体在移动，就得向Hittable类添加一个方法，让每个物体都注意到当前帧周期。这个方法将为所有动态物体提供一种在那一帧间设置它们运动的方式。

但我们接下来写的是一个简化模型。我们只会渲染单独一帧，隐式假设开始时间是0，结束时间是1。我们需要完成以下两个任务：

修改Camera类，支持随机生成时间在内的光线。
让球体类支持“移动”。

让摄像机模拟运动模糊

修改Camera类的get_ray()，让它返回时间在内的随机光线：

Ray get_ray(int i, int j) const
{
	Vec3 offset = sample_square();
	Point3 pixel_sample = pixel00_pos + ((i + offset.x()) * pixel_delta_u) + ((j + offset.y()) * pixel_delta_v);

	Point3 ray_origin = (defocus_angle <= 0) ? center : defocus_disk_sample();
	Vec3 ray_direction = pixel_sample - ray_origin;
	double ray_time = random_double();

	return Ray(ray_origin, ray_direction, ray_time);
}

添加移动的球体

接下来更新Sphere类，让球体动起来。需要定义两个球心，让球在time = 0时位于center1，在time = 1时位于center2，中间线性过渡：

class Sphere : public Hittable
{
public:
	// 静态球
	Sphere(const Point3& center, double radius, shared_ptr<Material> material)
		: center1(center), radius(std::fmax(0, radius)), material(material), is_moving(false)
	{}

	// 动态球
	Sphere(const Point3& center,  const Point3& center2, double radius, shared_ptr<Material> material)
		: center1(center), radius(std::fmax(0, radius)), material(material), is_moving(true)
	{
		center_vec = center2 - center1;
	}
	...
        
private:
	Point3 center1;
	double radius;
	shared_ptr<Material> material;
	bool is_moving;
	Vec3 center_vec;

	// 根据属于[0, 1]时间, 返回属于[center1, center2]线性插值的位置
	Point3 sphere_center(double time) const
	{
		return center1 + time * center_vec;
	}
};

另一种设计静态球的方法就是将所有球看成动态的，但静态球的起点和终点不变。接下来修改hit()方法，让其根据球体类型做改变：

bool hit(const Ray& r, Interval ray_t, HitRecord& rec) const override
{
	Point3 center = is_moving ? sphere_center(r.time()) : center1;
	// oc = 球心C - 光线原点Q
	Vec3 oc = center - r.origin();
    ...
}

支持求交时管理时间

别忘了更新所有材质的scatter()方法，让它们支持光线的时间变量：

scattered = Ray(rec.position, direction, r_in.time());

整合到一块

最后就可以进行整合了，我们直接用最终场景进行渲染，并将diffuse材质的球设置为动态球：

void book1Scene(HittableList& world, Camera& cam)
{
    auto ground_material = make_shared<Lambertian>(Color(0.5, 0.5, 0.5));
    world.add(make_shared<Sphere>(Point3(0, -1000, 0), 1000, ground_material));

    for (int a = -11; a < 11; a++) {
        for (int b = -11; b < 11; b++) {
            auto choose_mat = random_double();
            Point3 center(a + 0.9 * random_double(), 0.2, b + 0.9 * random_double());

            if ((center - Point3(4, 0.2, 0)).length() > 0.9) {
                shared_ptr<Material> sphere_material;

                if (choose_mat < 0.8) {
                    // diffuse
                    auto albedo = Color::random() * Color::random();
                    sphere_material = make_shared<Lambertian>(albedo);
                    auto center2 = center + Vec3(0, random_double(0, 0.5), 0);
                    world.add(make_shared<Sphere>(center, center2, 0.2, sphere_material));
                }
                ...
    }
	...

    cam.aspectRadio = 16.0 / 9.0;
    cam.imgWidth = 400;
    cam.samples_per_pixel = 100;
    cam.max_depth = 50;

    cam.vfov = 20;
    cam.lookFrom = Point3(13, 2, 3);
    cam.lookAt = Point3(0, 0, 0);
    cam.vup = Vec3(0, 1, 0);

    cam.defocus_angle = 0.6;
    cam.focus_dist = 10.0;
}

最终成果如下：

层次包围体结构BVH

光线-物体求交是光线追踪器的性能瓶颈，并且运行时间随着物体数量线性增加。但这个过程是在相同场景的重复搜索，所以我们让它成为基于二叉树的对数搜索。由于我们朝相同场景中发射成百上千万条光线，可以排序场景中的对象，然后每次光线求交都是一个次线性搜索。可以通过划分空间或划分对象的方式实现，后者相对容易一些。

核心思想

为一堆图元创建包围盒的核心思想就是找到一个能完全包裹住所有对象的包围体。例如有个包住10个对象的球体，任何不与这个球相交的光线都不会和这10个对象相交，因此代码的结构如下：

if (光线打到包围盒)
    return 是否打到被包围的物体
else
    return false

在这里我们不划分场景空间，而是为每一个对象准备一个包围盒，即使它们会重合。

包围体的层次结构

为了让搜索变成次线性的，我们需要构造BVH。举个例子，如果我们将物体分为红蓝两组，且用矩形包围盒，会有如下结构：

注意蓝色和红色包围盒被紫色的包围盒包含，而且它们可能重叠&没有规律。因此右边树的左右子节点也没有规律，代码结构如下：

if (光线打到紫色包围盒)
{
    hit0 = 光线打到蓝色包围盒的物体;
    hit1 = 光线打到红色包围盒的物体;
    if (hit0 || hit1)
        return true 和 打到物体的详细信息
}
return false;

轴对齐包围盒AABB

为了得到高效的BVH，需要高效的划分方法和光线与包围盒的相交方法。光线与包围盒的相交应该很快，并且包围盒需要紧密包含对象。对于大多数模型使用轴对齐包围盒AABB就够了。

以二维的AABB为例，它由x和y两段区间相交而成：

接下来看看如何计算光线是否与AABB相交。首先是x轴，假设光线和x轴区间相交于和两点（如果光线和x轴平行，这两点未定义）：

回想光线的数学定义：将其拆分为x，y，z三个分量，然后讨论x轴的情况。这里让，有：通过移项便能求出：同理求出：将一维转换到二维，三维的关键点就是：如果光线和包围盒的所有平面相交，那么这几对的区间会重叠。例如下图的2D情况，蓝色和绿色区间重叠的前提条件是光线和2D包围盒相交：

光线和AABB的求交

光线和二维AABB是否相交的代码如下：

tx区间 = 在x轴计算相交(ray, x0, x1);
ty区间 = 在y轴计算相交(ray, y0, y1);
return 区间是否重叠(tx区间, ty区间);

这看起来挺简单，且容易拓展到三维情况：

tx区间 = 在x轴计算相交(ray, x0, x1);
ty区间 = 在y轴计算相交(ray, y0, y1);
tz区间 = 在z轴计算相交(ray, z0, z1);
return 区间是否重叠(tx区间, ty区间, tz区间);

回想一下1D情况下的和：为了避免区间边界问题，建议写成如下形式：但仍无法避免等问题，我们将在以后讨论。

接下来看看上面代码的区间是否重叠()怎么写：

bool overlaps(t_interval1, t_interval2)
{
    t_min = max(t_interval1.min, t_interval2.min);
    t_max = min(t_interval1.max, t_interval2.max);
    return t_min < t_max;
}

如果在计算过程中出现NaN，结果会返回false，所以我们需要处理好边界条件。为了处理好边界条件，先在Interval类中添加方法expand(delta)，该方法将为整个区间均匀扩大delta：

// 将整个区间均匀扩大delta
Interval expand(double delta) const
{
	double padding = delta / 2;
	return Interval(min - padding, max + padding);
}

然后实现包围盒类AABB：

#pragma once

#include "../rtweekend.h"

class AABB
{
public:
	Interval x, y, z;

	// 创建一个空AABB
	AABB() = default;
	// 创建一个区间分量为x, y, z的AABB
	AABB(const Interval& x, const Interval& y, const Interval& z)
		: x(x), y(y), z(z) {}
	// 用两个点创建一个AABB
	AABB(const Point3& a, const Point3& b)
	{
		x = (a[0] <= b[0]) ? Interval(a[0], b[0]) : Interval(b[0], a[0]);
		y = (a[1] <= b[1]) ? Interval(a[1], b[1]) : Interval(b[1], a[1]);
		z = (a[2] <= b[2]) ? Interval(a[2], b[2]) : Interval(b[2], a[2]);
	}

	// 根据索引返回对应坐标轴区间, 0-x, 1-y, 2-z
	const Interval& axis_interval(int n) const
	{
		if (n == 1)	return y;
		if (n == 2)	return z;
		return x;
	}

	// 判断光线和AABB求交
	bool hit(const Ray& r, Interval ray_t) const
	{
		const Point3& ray_orig = r.origin();
		const Vec3& ray_dir = r.direction();

		for (int axis = 0; axis < 3; ++axis)
		{
			const Interval& ax = axis_interval(axis);
			const double adinv = 1.0 / ray_dir[axis];

			double t0 = (ax.min - ray_orig[axis]) * adinv;
			double t1 = (ax.max - ray_orig[axis]) * adinv;

			if (t0 < t1)
			{
				if (t0 > ray_t.min)	ray_t.min = t0;
				if (t1 < ray_t.max)	ray_t.max = t1;
			}
			else
			{
				if (t1 > ray_t.min)	ray_t.min = t1;
				if (t0 < ray_t.max)	ray_t.max = t0;
			}

			if (ray_t.max <= ray_t.min)
			{
				return false;
			}
		}

		return true;
	}
};

为物体添加包围盒

现在需要为所有种类物体添加计算包围盒的方法，然后将为这些图元建立BVH。对于球体这种独立的图元，将作为BVH的叶子节点存在。

默认的AABB为空，因此一些物体可能有空的包围盒。例如，考虑一个没有子节点的HittableList。不过不用担心计算问题，因为Interval类会解决。最后回想一下有的物体可能是运动的，那它的包围盒就得是整个运动范围，从time = 0到time = 1。

首先，先给Hittable类添加获取AABB的虚函数：

class Hittable 
{
public:
	virtual ~Hittable() = default;

	virtual bool hit(const Ray& r, Interval ray_t, HitRecord& rec) const = 0;

	virtual AABB bounding_box() const = 0;
};

然后是静态球的AABB：

class Sphere : public Hittable
{
public:
	// 静态球
	Sphere(const Point3& center, double radius, shared_ptr<Material> material)
		: center1(center), radius(std::fmax(0, radius)), material(material), is_moving(false)
	{
		Vec3 rvec = Vec3(radius, radius, radius);
		bbox = AABB(center1 - rvec, center1 + rvec);
	}

	...
        
	AABB bounding_box() const override { return bbox; }

private:
	Point3 center1;
	double radius;
	shared_ptr<Material> material;
	bool is_moving;
	Vec3 center_vec;
	AABB bbox;

	...
};

对于动态球，我们想要它的整个运动范围。可以分别获取它运动开始和结束时的两个AABB，然后合并两个AABB：

// 动态球
Sphere(const Point3& center,  const Point3& center2, double radius, shared_ptr<Material> material)
	: center1(center), radius(std::fmax(0, radius)), material(material), is_moving(true)
{
	Vec3 rvec = Vec3(radius, radius, radius);
	AABB box0(center1 - rvec, center1 + rvec);
	AABB box1(center2 - rvec, center2 + rvec);
	bbox = AABB(box0, box1);

	center_vec = center2 - center1;
}

然后我们需要实现AABB新的构造函数，在此要先添加一个Interval类的构造函数：

// 合并两个区间
Interval(const Interval& a, const Interval& b)
{
	min = a.min <= b.min ? a.min : b.min;
	max = a.max >= b.max ? a.max : b.max;
}

接下来就能实现新的AABB构造函数了：

// 合并两个包围盒
AABB(const AABB& box0, const AABB& box1)
{
	x = Interval(box0.x, box1.x);
	y = Interval(box0.y, box1.y);
	z = Interval(box0.z, box1.z);
}

对于三角形组成的模型，首先在Model类中添加AABB私有变量和构建它的build()方法：

class Model : public Hittable
{
public:
	Model(const std::vector<Triangle>& triangles, shared_ptr<Material> material)
		: triangles(triangles), material(material)
	{
		build_aabb();
	}
	Model(const std::filesystem::path& filename, shared_ptr<Material> material)
		: material(material)
	{
		...
            
		build_aabb();
	}

    ...

	AABB bounding_box() const override { return bbox; }

private:
	AABB bbox;
	std::vector<Triangle> triangles;
	shared_ptr<Material> material;

	// 构建包围盒
	void build_aabb()
	{
		for (const auto& triangle : triangles)
		{
			bbox = AABB(bbox, triangle.bounding_box());
		}
	}
};

然后修改Triangle类，添加包围盒和相关方法：

class Triangle : public Hittable
{
public:
	Triangle(const Point3& p0, const Point3& p1, const Point3& p2,
			 const Vec3& n0, const Vec3& n1, const Vec3& n2,
			 shared_ptr<Material> material)
		: p0(p0), p1(p1), p2(p2), n0(n0), n1(n1), n2(n2), material(material)
	{
		build_aabb();
	}

	Triangle(const Point3& p0, const Point3& p1, const Point3& p2, 
			 shared_ptr<Material> material)
		: p0(p0), p1(p1), p2(p2), material(material) 
	{
		...
            
		build_aabb();
	}

	...

	AABB bounding_box() const override { return bbox; }

private:
	Point3 p0, p1, p2;	// 三角形的顶点位置信息
	Vec3 n0, n1, n2;	// 三角形的顶点法线信息
	shared_ptr<Material> material;
	AABB bbox;

	void build_aabb()
	{
		bbox = AABB(p0, p1, p2);
	}
};

最后添加AABB类的构造方法即可：

// 用三个点创建一个AABB
AABB(const Point3& a, const Point3& b, const Point3& c)
{
	auto [minX, maxX] = std::minmax({a[0], b[0], c[0]});
	x = Interval(minX, maxX);
	auto [minY, maxY] = std::minmax({a[1], b[1], c[1]});
	y = Interval(minY, maxY);
	auto [minZ, maxZ] = std::minmax({a[2], b[2], c[2]});
	z = Interval(minZ, maxZ);
}

为HittableList添加包围盒

别忘了给HittableList类添加包围盒，它的包围盒将随物体列表的增加而扩大：

public:
    void add(shared_ptr<Hittable> object) 
    {
        objects.push_back(object);
        bbox = AABB(bbox, object->bounding_box());
    }
...
	AABB bounding_box() const override { return bbox; }

private:
	AABB bbox;

BVH节点类

BVH也属于Hittable，它是一个容器，但可以去询问光线是否击中它。BVH有两种设计方式，第一种是设计节点类和树类，第二种是二者结合。这里选择第二种，对于hit()方法，只需检查该节点的包围盒是否被击中，是的话就检查子节点的包围盒。

BVHNode类设计如下：

#pragma once

#include "../rtweekend.h"
#include "aabb.hpp"
#include "../ray/hittable.hpp"
#include "../ray/hittableList.hpp"

class BVHNode : public Hittable
{
public:
	BVHNode(HittableList list) : BVHNode(list.objects, 0, list.objects.size()) {}

	BVHNode(std::vector<shared_ptr<Hittable>>& objects, size_t start, size_t end)
	{
		// TODO
	}

	bool hit(const Ray& r, Interval ray_t, HitRecord& rec) const override
	{
		if (!bbox.hit(r, ray_t))
		{
			return false;
		}

		bool hit_left = left->hit(r, ray_t, rec);
		bool hit_right = right->hit(r, Interval(ray_t.min, hit_left ? rec.t : ray_t.max), rec);

		return hit_left || hit_right;
	}

	AABB bounding_box() const override { return bbox; }

private:
	shared_ptr<Hittable> left;
	shared_ptr<Hittable> right;
	AABB bbox;
};

建立BVH

对于任何效率方面的数据结构，最复杂的就是建立它。我们在构造函数里建立BVH。无论列表中的对象如何划分成两个子列表，hit()都能正常工作。如果划分的好，效率会很高。这里将按轴的长度划分：

随机选一个轴
用std::sort()排列图元
每个子树放一半图元

当传入的列表只有两个元素，会在每个子树各放一个元素然后结束递归。遍历BVH的算法应该平滑且不必检查空指针，因此只有一个元素时会在两个子树都添加该元素。建立BVH的代码如下：

BVHNode(std::vector<shared_ptr<Hittable>>& objects, size_t start, size_t end)
{
	int axis = random_int(0, 2);
	auto comparator = (axis == 0) ? box_x_compare
					: (axis == 1) ? box_y_compare
								  : box_z_compare;

	size_t object_span = end - start;
	if (object_span == 1)
	{
		left = right = objects[start];
	}
	else if (object_span == 2)
	{
		left = objects[start];
		right = objects[start + 1];
	}
	else
	{
		std::sort(std::begin(objects) + start, std::begin(objects) + end, comparator);
		
		size_t mid = start + object_span / 2;
		left = make_shared<BVHNode>(objects, start, mid);
		right = make_shared<BVHNode>(objects, mid, end);
	}

	bbox = AABB(left->bounding_box(), right->bounding_box());
}

接下来在rtweekend.h编写生成随机整数的random_int()：

inline int random_int(int min, int max)
{
	return static_cast<int>(random_double(min, max + 1));
}

对于无限大平面这类的图元不需要包围盒，如果要添加这个东西，需要再改一下构造函数。

包围盒比较函数

接下来需要实现给std::sort()用的包围盒比较函数。首先我们创建一个通用的比较器如下：

static bool box_compare(const shared_ptr<Hittable> a, const shared_ptr<Hittable> b, int axis_idx)
{
	Interval a_axis_interval = a->bounding_box().axis_interval(axis_idx);
	Interval b_axis_interval = b->bounding_box().axis_interval(axis_idx);
	return a_axis_interval.min < b_axis_interval.min;
}

然后用这个通用的比较器实现三个轴的专用比较器：

static bool box_x_compare(const shared_ptr<Hittable> a, const shared_ptr<Hittable> b)
{
    return box_compare(a, b, 0);
}

static bool box_y_compare(const shared_ptr<Hittable> a, const shared_ptr<Hittable> b)
{
    return box_compare(a, b, 1);
}

static bool box_z_compare(const shared_ptr<Hittable> a, const shared_ptr<Hittable> b)
{
    return box_compare(a, b, 2);
}

接下来就能用BVH了：

world = HittableList(make_shared<BVHNode>(world));

一些优化

在BVH的构造函数中，选择随机一个轴进行划分，实际上可以选最长轴进行划分。

首先，先为BVH创建一个空AABB，然后拓展其大小直至能够容纳所有对象。然后就能找到AABB的最长轴，然后用最长轴进行划分了。修改BVH构造函数如下：

// 构造容纳所有对象的AABB
bbox = AABB::empty;
for (size_t idx = start; idx < end; ++idx)
{
	bbox = AABB(bbox, objects[idx]->bounding_box());
}

// 按最长轴划分
int axis = bbox.longest_axis();
auto comparator = (axis == 0) ? box_x_compare : (axis == 1) ? box_y_compare : box_z_compare;

接下来完善AABB类，添加相关方法：

class AABB
{
    ...
        
	// 返回最长轴
	int longest_axis() const
	{
		if (x.size() > y.size())
		{
			return x.size() > z.size() ? 0 : 2;
		}
		else
		{
			return y.size() > z.size() ? 1 : 2;
		}
	}

	static const AABB empty, universe;
};

const AABB AABB::empty = AABB(Interval::empty, Interval::empty, Interval::empty);
const AABB AABB::universe = AABB(Interval::universe, Interval::universe, Interval::universe);

这样就会渲染地稍微快一点了。

参考资料

Ray Tracing: The Next Week