从0到实时光线追踪(一):基本数据结构

这一系列教程更倾向于从工程领域展开思路,换句话说就是会更多地进行代码上的实现。同时限制于笔者自身水平,在这里我姑且假设读者都有基本的现代图形学知识,换句话说就是已经完成了GAMES101的学习,以至于不必过多地在最基础的部分上做过多讲解,拉长篇幅。

本文的目的是构造一个轻量级的实时光线追踪渲染器。在本章,我们会构造渲染器的基本数据结构,以支撑随后的整体渲染架构。

我们的架构基于python的Taichi框架,因此或许在开始之前,你需要掌握一定的python语法,安装3.8及以上版本的python解释器,并:

1
npm pip install taichi

在随后的编写过程中,我们默认按这样的方式调用了Taichi框架:

1
2
import taichi as ti
import taichi.math as tm

Ray

在实时光线追踪中,我们通常会使用光栅化-光线追踪同时进行的混合渲染方法进行渲染,这种渲染管线也被称为混合渲染管线(hybrid renderer pipeline)。因此,我们需要实现光线追踪中的光线类。在现阶段,我们只需要实现一个具有原点和方向的简单射线就好:

其中表示物体的原点,表示物体的方向。同时,由于光线需要表达其传播的色彩,我们也需要记录射线的颜色。现在,我们暂时使用三维向量储存颜色,而在之后我们将实现自己的颜色类与混色算法。

1
2
3
4
5
6
7
8
9
10
@ti.dataclass
class Ray:
origin: vec3
direction: vec3
color: vec3

def __init__(self, ori, dire, color):
self.origin = ori
self.direction = dire
self.color = color

Transform

基本定义

在一般的图形学渲染器中,我们使用Transform来记录物体的位置、旋转与尺寸,而在对一个物体进行位移、旋转和缩放物体时,只需要调整其Transform就好。通常,我们会首先使用三维的列向量形式记录位移、旋转与缩放的值:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
from taichi.math import vec2,vec3,vec4,mat3,mat4

@ti.dataclass
class Transform:
# data
position: vec3
rotation: vec3
scale: vec3

def __init__(self, p, r, s):
self.position = p
self.rotation = r
self.scale = s

在这之中,物体的位置和缩放很容易处理,而旋转则似乎难以下手。事实上,我们可以将旋转角度转写为矩阵的形式,这样就可以通过矩阵与列向量之间的乘法来计算得到物体旋转后的位置信息了。

这里并不会阐述具体的推导过程。不过,只要据此进行计算,就能得到旋转角度对应的旋转矩阵。再将旋转矩阵与物体位置信息相乘,就得到了物体旋转后的新位置信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def update_rotate(p, self) -> mat3:
s = sin(self.rotation)
c = cos(self.rotation)
rotation_matrix = mat3(vec3(c.z, s.z, 0),
vec3(-s.z, c.z, 0),
vec3(0, 0, 1)) @ \
mat3(vec3(c.y, 0, -s.y),
vec3(0, 1, 0),
vec3(s.y, 0, c.y)) @ \
mat3(vec3(1, 0, 0),
vec3(0, c.x, s.x),
vec3(0, -s.x, c.x))
return rotation_matrix
# 在Taichi中,我们使用A @ B表示A和B两个矩阵之间的矩阵乘法
# 如果使用A * B,则代表A和B矩阵各个数位上的数分别相乘所得到的新矩阵

显而易见的是,这样的矩阵运算会耗费大量的时间,因此我们通常会在Transform中同时存储物体的旋转矩阵,以达到节约时间开销的目的。

1
2
3
4
5
def __init__(self, p, r, s):
self.position = p
self.rotation = r
self.scale = s
self.rotation_matrix = self.rotate()

此时,我们的Transform即为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
class Transform:
# 基本定义
position: vec3
rotation: vec3
scale: vec3
rotation_matrix: mat3

def __init__(self, p, r, s):
self.position = p
self.rotation = r
self.scale = s
self.rotation_matrix = self.rotate()

def update_rotate(self) -> mat3:
s = sin(self.rotation)
c = cos(self.rotation)
rotation_matrix = \
mat3(vec3(c.z, s.z, 0),
vec3(-s.z, c.z, 0),
vec3(0, 0, 1)) @ \
mat3(vec3(c.y, 0, -s.y),
vec3(0, 1, 0),
vec3(s.y, 0, c.y)) @ \
mat3(vec3(1, 0, 0),
vec3(0, c.x, s.x),
vec3(0, -s.x, c.x))
return rotation_matrix

使用四元数类

然而,使用欧拉角完成的旋转往往会带来两大问题——一是欧拉死锁,二是旋转插值无法进行。欧拉死锁指在绕一定方向进行一定的欧拉旋转后,物体的某一轴向将会与另一轴向相重合,导致物体的旋转维度降低,旋转插值无法进行则是因为欧拉角旋转并没有指定方向:举个例子,如果我们先将一个物体顺时针旋转了359度,再将其逆时针旋转了361度,那在进行第二步旋转时,究竟是应该将物体逆时针旋转361度,还是应该将其顺时针旋转2度呢?

虽然在离线光线追踪渲染中,这一问题十分容易避免——使用运动模糊与大批量的单帧渲染掩盖这一问题——但是在实时的渲染中,我们需要使用更恰当的方法解决它,而这一方法就是大名鼎鼎的四元数旋转(Quaternion)。

四元数可以看作是复数的拓展,它由一个实部和三个虚部构成。一般而言写作:

其中,为实数部分,为基本虚数单位。四元数旋转的核心思想是将旋转操作转换为四元数乘法,并通过将四元数乘法转换为矩阵运算的形式来实现旋转操作。我们通常通过这样的步骤进行四元数旋转:

  • 将旋转角度和旋转轴转换为四元数表示,其中旋转角度为弧度制,旋转轴必须是单位向量
  • 计算初始四元数,将旋转轴乘以得到虚数部分,将本身作为实数部分,得到形如上文形式的四元数
  • 将需要旋转的向量转换为四元数表示,其中实数部分为0,虚数部分使用分量表示。
  • 计算旋转后的四元数,即将初始四元数和向量四元数相乘。
  • 将旋转后的四元数转回欧拉表示。

通常,我们使用一个向量和一个变量来表示四元数:

1
2
3
4
5
6
7
8
@ti.dataclass
class Quaternion:
v: vec3
w: float

def __init__(self, v: vec3, w: float):
self.v = v
self.w = w

同时,四元数的四则运算遵循这样的法则:

  • 对于四元数之间的加减法,直接将对应的实部与虚部相加减即可;
  • 四元数与实数相乘,直接将四元数的各部与实数分别相乘即可;
  • 四元数与四元数相乘,则有:

因此,我们可以重载四元数的四则运算:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def __add__(self, other):
return Quaternion(self.v + other.v, self.w + other.w)

def __sub__(self, other):
return Quaternion(self.v - other.v, self.w - other.w)

def __mul__(self, other):
if isinstance(other, Quaternion):
return Quaternion(
cross(self.v, other.v) + other.v * self.w + self.v * other.w,
self.w * other.w - dot(self.v, other.v)
)
elif isinstance(other, (int, float)):
return Quaternion(self.v * other, self.w * other)

同时,在进行旋转前,我们需要将四元数转化为单位四元数:

1
2
3
4
def normalize(self):
p = self.v.x ** 2 + self.v.y ** 2 + self.v.z ** 2 + self.w ** 2
inv_p = 1 / sqrt(p)
return Quaternion(self.v * inv_p, self.w * inv_p)

统合以上的所有代码,我们就得到了四元数类:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
from taichi.math import vec3, vec4
from taichi.math import cross, dot, sqrt


@ti.dataclass
class Quaternion:
v: vec3
w: float

def __init__(self, v: vec3, w: float):
self.v = v
self.w = w

# 运算符重载
def __add__(self, other):
return Quaternion(self.v + other.v, self.w + other.w)

def __sub__(self, other):
return Quaternion(self.v - other.v, self.w - other.w)

def __mul__(self, other):
if isinstance(other, Quaternion):
return Quaternion(
cross(self.v, other.v) + other.v * self.w + self.v * other.w,
self.w * other.w - dot(self.v, other.v)
)
elif isinstance(other, (int, float)):
return Quaternion(self.v * other, self.w * other)

def __str__(self):
return "[{}+{}i+{}j+{}k]".format(self.w, self.v.x, self.v.y, self.v.z)

# 方法
def normalize(self):
p = self.v.x ** 2 + self.v.y ** 2 + self.v.z ** 2 + self.w ** 2
inv_p = 1 / sqrt(p)
return Quaternion(self.v * inv_p, self.w * inv_p)

不过,由于目前还没有涉及到动画的问题,我们暂时不会应用四元数类,因为在静止不动的单帧渲染中欧拉角是一个十分方便且快捷的工具。当我们进行场景动画处理时,四元数会变得十分常用,因此提前对它有一个了解是很重要的一件事。

Material

除了位置信息,一个物体还需要具有材质信息。我们使用一个最简单的真实感物理材质(PBR)来完成我们的初版渲染器,这个材质包含以下信息:

  • Albedo:固有色,反映了物体本身的颜色

  • Roughness:粗糙度,反映了物体表面的粗糙程度

  • Metallic:金属度,反映了物体表面的光泽程度

  • Transmission:透明度,反映了物体是否透明

  • IOR:折射率,在物体发生折射时,表示折射角的方向

  • Emission:自发光,由两部分组成,一是自发光强度,二是自发光颜色。

在现阶段,我们只需要给予材质一个基本的定义即可。在随后的渲染过程中,我们才会对这些参数进行应用。基于此,我们可以写出一个简单的类:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
@ti.dataclass
class Material:
albedo: vec3
roughness: float
metallic: float
transmission: float
IOR: float
emission_intensity: float
emission_color: vec3

def __init__(self, albedo, roughness, metallic, transmission=0, IOR=1,
emission_intensity=0, emission_color=vec3(0)):
self.albedo = albedo
self.roughness = roughness
self.metallic = metallic
self.transmission = transmission
self.IOR = IOR
self.emission_intensity = emission_intensity
self.emission_color = emission_color

BoundingBox

一般来说,直接处理物体与物体、物体与光线等之间的碰撞问题的开销十分巨大——这要求我们逐个顶点地检测是否发生了碰撞,而大多物体都有成千上万个顶点。为了减少程序的计算量,提升渲染速度,我们通常会使用到包围盒(BoundingBox)。包围盒是一个封装了物体的最小立方体、最小矩形或最小球体。在检测物体碰撞时,如果两个物体的包围盒都没有碰撞,那我们就可以认为这两个物体本身没有发生碰撞,通过检测包围盒,就可以避免大量的无效计算了。

常见的包围盒有:

  • AABB包围盒(Axis-aligned bounding box):是一种与坐标轴平行的矩形框,最简单、最常见的包围盒类型。
  • OBB包围盒(Oriented bounding box):是一个可旋转的长方体,其面不一定与坐标轴平行。由于 OBB 可以更好地适应物体的外形,因此它能提供更真实的碰撞检测结果。但同时也需要更高的计算成本。
  • 包围球(Bounding Sphere):一个包围物体的最小球体,检测物体之间的相对距离十分方便。

在目前阶段,我们仅使用一个AABB包围盒。对于AABB包围盒,我们只需要定义两个顶点即可。

1
2
3
4
5
6
7
8
9
@ti.dataclass
class BoundingBox:
p_min: vec3
p_max: vec3

def __init__(self, p_min, p_max):
self.p_min = p_min
self.p_max = p_max

包围盒需要定义几个简单的函数。首先是,我们需要将一个点加入已有的包围盒,这只需要比较新加入的点坐标与包围盒边界点坐标的大小即可。Taichi为我们提供了十分方便的方式来更新包围盒的边界点:

1
2
3
def add(self, p: vec3):
self.p_min = min(self.p_min, p)
self.p_max = max(self.p_max, p)

同时,有时我们也需要合并两个包围盒。这同样只需要比较两个包围盒各自的最小坐标与最大坐标即可。美观起见,我们将方法定义在类外。

1
2
def bounding_box_merge(a, b):
return BoundingBox(min(a.p_min, b.p_min), max(a.p_max, b.p_max))

除了修改包围盒的范围之外,我们也需要判断一个点是否在包围盒的范围内。令人兴奋的是,Taichi提供了甜美的语法糖来帮助我们完成这个简洁却十分优雅的方法:

1
2
3
4
5
def contain(self, p: vec3) -> bool:
s = (self.p_min < p < self.p_max)
return sum(s) == 3
# 在Taichi中,当比较两个n维Vector的大小时,将会得到一个长度为n的列表,记录了在各个维度比较大小的结果
# 显而易见的是,如果在3个维度比较大小的结果都是True,那么点就在包围盒中

最后,我们需要实现判断射线是否能与包围盒相交。判断包围盒与射线的相交十分简单,我们只需要知道射线进入包围盒的时间和离开包围盒的时间,若满足:

即可。射线的进入时间就是射线分别进入三轴的最大时间,离开时间就是射线分别退出三轴的最小时间:

据此我们就可以得到检测射线与包围盒相交的函数。同样,基于Taichi提供的语法糖,这一函数的实现可以变得极其优雅:

1
2
3
4
5
6
7
8
9
10
11
12
@ti.func
def ray_inter_bouding_box(r: Ray, b: BoundingBox) -> bool:
# o+dt=s --> t=(s-o)/d
t_a = (b.p_min - r.origin) / (r.direction + vec3(1e-5)) # avert the zero-division
t_b = (b.p_min - r.origin) / (r.direction + vec3(1e-5))
t_enter_xyz = min(t_a, t_b)
t_exit_xyz = max(t_a, t_b)

t_enter = max(t_enter_xyz.x, t_enter_xyz.y, t_enter_xyz.z)
t_exit = min(t_exit_xyz.x, t_exit_xyz.y, t_exit_xyz.z)

return t_enter < t_exit and t_exit_xyz >= 0

SDFObject

有了Transform和Material,自然需要有与之配套的物体。我们首选使用有向距离场物体进行场景测试。有向距离场物体(SDF Object)是一类使用公式表示的隐式几何物体,将任意一个点带入一个有向距离场公式,若得到的结果大于0,说明点p在这个物体的外面;若得到的结果等于0,说明点p在物体表面上;若得到的结果小于0,说明点p在物体的内部。

显然,有多少种公式,就能得到多少种不一样的有向距离场形状物体,因此我们的SDFObject类中首先要定义物体的类型type。随后,我们需要给物体绑定一个Transform,表示物体的位置信息。最后,物体也需要带有一个包围盒。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
SPHERE = 0
CUBE = 1


@ti.dataclass
class SDFObject:
type: int
trans: Transform
mtl: Material
bound: BoundingBox

def __init__(self, type, trans, mtl):
self.type = type
self.trans = trans
self.mtl = mtl

作为测试对象,我们给出CUBE和SPHERE的有向距离场公式,并使用一个函数统一调度两个公式:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
@ti.func
def sdf_sphere(p: tm.vec3, r: float) -> float:
return tm.length(p) - r


@ti.func
def sdf_cube(p: tm.vec3, b: tm.vec3) -> float:
q = abs(p) - b
return tm.length(tm.max(q, 0)) + tm.min(tm.max(q.x, q.y, q.z), 0)

@ti.func
def sdf_dis(p: tm.vec3, obj: SDFObject) -> float:
dis = 0
if obj.type == CUBE:
dis = sdf_cube(p, obj.trans.scale)
elif obj.type == SPHERE:
dis = sdf_sphere(p, obj.trans.scale)

return dis

同时,当物体的Transform发生改变时,需要同时改变其包围盒,因此我们需要实现一个简单的包围盒更新方法。在通过坐标更新包围盒时,为了修正精度损失带来的包围盒大小误差,我们将略微对包围盒进行扩大。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def update_bound(self):
pos = self.trans.position
scale = self.trans.scale
rot_mat = self.trans.rotation_matrix
if self.type == SPHERE:
self.bound.p_min = vec3(pos - scale.x)
self.bound.p_max = vec3(pos + scale.x)
elif self.type == CUBE:
a = ti.abs(scale)
b = -ti.abs(scale)
p_min = vec3(1e9, 1e9, 1e9)
p_max = vec3(-1e9, -1e9, -1e9)
for i in range(8):
v = vec3(
a.x if (i & 1) == 0 else b.x,
a.y if (i & 2) == 0 else b.y,
a.z if (i & 4) == 0 else b.z
)
v = rot_mat @ v
p_min = min(p_min, v)
p_min = min(p_min, v)
p_max = max(p_max, v)

self.bound.p_min = p_min * 1.15 + self.trans.position
self.bound.p_max = p_max * 1.15 + self.trans.position

最后,为了美观起见,我们将所有的数据更新方法封装为一个方法,并在每次模型的数据更新后调用它。

1
2
3
def update_state(self):
self.trans.update_rotate()# 更新旋转矩阵
self.update_bound()# 更新包围盒

有关更一般的.obj格式模型与.fbx格式模型等的设计将会在完成渲染器的光栅化处理后给出。

Camera

我们需要使用一个相机来渲染场景。虽然Taichi已经封装好了一个可用的相机,但为了加强可自定义程度,并为随后实现光场相机做准备,在这里我们仍然选择自定义一个相机。

通常情况下,一个相机都含有三个基本参数:position用于记录相机的位置,lookat用于记录相机的朝向,或者说相机当前指向的点,up指定相机的上方。同时,相机的aspect指定了相机平面的宽高比,fov制定了相机的视野角度。以上的五个参数共同构成了一个在其它渲染器中会提供的基本相机。而为了实现具备物理特征的相机,我们会额外定义光圈大小aperture用于控制景深,focus用于控制对焦距离。

1
2
3
4
5
6
7
8
9
@ti.dataclass
class Camera:
position: vec3
lookat: vec3
up: vec3
fov: float
aspect: float
aperture: float # 光圈大小
focus: float # 焦距

在场景中,通常只会存在一个相机。因此我们可以选择在定义 Camera 类之后直接生成这个相机。同时,为了更好地管理相机的属性,我们可以新建一个配置文件 RenderConfig.py,从中读取相机的各个属性。

当我们需要对一个像素进行渲染时,我们首先将该像素在屏幕上的位置传递给相机。通过相机的视角大小和宽高比,我们可以计算出相机平面的大小,并生成一条从相机位置朝向前方的射线。这条射线在光线追踪和光栅化中都起到重要作用。因此,我们需要仔细地计算相机平面上任意一点所对应的射线方向。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
@ti.func
def get_direct_ray(p: vec3, l: vec3, u: vec3, uv: vec2, color: vec3) -> Ray:
position = p
lookat = l
up = u

# compute the half angle of the FOV of the camera
theta = radians(scene_camera.fov)
half_height = tan(theta * 0.5)
half_width = scene_camera.aspect * half_height

# compute the camera coordinate system
z = normalize(position - lookat)
x = normalize(cross(up, z))
y = cross(z, x)

# compute the start position of the ray
start_pos = position - half_width * x - half_height * y - z

# compute the horizontal and vertical size of the plane
horizontal = 2.0 * half_width * x
vertical = 2.0 * half_height * y

# get ray
ro = position
rp = start_pos + uv.x * horizontal + uv.y * vertical
rd = normalize(rp - ro)

return Ray(ro, rd, color)

然而,这一射线方向虽说能真实地向场景中发射光线,但是这一光线发射方式却没有考虑焦距、光圈等物理参量对光线的影响,因此并不能用于物理相机。在这里,我们提供一个额外的光线传播函数,用以模拟物理相机的各种参量。我们采用Thin Lens模型模拟物理相机的焦距与光圈信息。

我们首先根据相机的光圈aperture计算出相机的光圈半径:

随后使用随机向量和光圈半径计算出光圈上任意一点的位置,使得采样到的位置满足径向分布函数。在计算完光圈采样点的偏离量之后,就可以借助光源与相机光轴的焦点作为射线的起点,从而得到每条光线从相机出发的方向向量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
@ti.func
def get_ray(p: vec3, l: vec3, u: vec3, uv: vec2, color: vec3):
position = p
lookat = l
up = u
focus = scene_camera.focus

# compute the half angle of the FOV of the camera
theta = radians(scene_camera.fov)
half_height = tan(theta * 0.5)
half_width = scene_camera.aspect * half_height

# compute the camera coordinate system
z = normalize(position - lookat)
x = normalize(cross(up, z))
y = cross(z, x)

# compute the start position of the ray
start_pos = position - half_width * x * focus - half_height * y * focus - z * focus

# compute the horizontal and vertical size of the plane
horizontal = 2.0 * half_width * x * focus
vertical = 2.0 * half_height * y * focus

# compute the ray direction from the camera to the target pixel
lens_radius = scene_camera.aperture * 0.5
a = ti.random()
unit_disk_random = sqrt(ti.random()) * vec2(sin(a), cos(a))
rud = lens_radius * unit_disk_random
offset = x * rud.x + y * rud.y

# get ray
ro = position + offset
rp = start_pos + uv.x * horizontal + uv.y * vertical
rd = normalize(rp - ro)

return Ray(ro, rd, color)

这就是所有的基本类了。其它诸如之类的大多数数学容器都已经在Taichi中有所实现。我们的第一个小目标是渲染出光栅化的静态CornellBox场景。

到目前为止,我们已经完成了以下工作:

[i835xy.png]