引子

最近在照着 Ray Tracing In One Weekend 系列用 Python 实现一个蒙特卡洛光追器（代码在这里）。第一本书的内容已经写完，但是由于 Python 感人的性能，我决定在学第二本书之前先把这个尚不完备的光追器进行一些性能优化。终于，优化后的性能可以让我在能够接受的时间内渲染出一张足够炫的图了。

然而，在渲染第一本书最后一章的场景时，我碰到了一个神秘的 Bug。运行的时候，虽然不会出现 Exception 停止渲染，但是 Numpy 会出现这样的运行时警报：

1...\InOneWeekend_optimized\utils\vec3.py:95: RuntimeWarning: invalid value encountered in sqrt
2  normal * (-np.sqrt(1 - r_out_parallel.length_squared()))

这里其实一共会出现 8 个 RuntimeWarning，我只附上了引起连锁问题的第一个警报。

这是一个相当偶发的 Bug，在一整张 1080P 分辨率、大约几十的采样数和递归深度的图渲染的过程中，这个问题只会出现一到两次。而发生问题的这个位置是在 Vec3 类中计算透明材质折射的位置出现的，这个函数的调用量绝对巨大。

我其实完全可以忽略这个问题，因为实在偶发，并且也不会影响结果。但是去他的，我花了整整两天时间找到了问题的根源。

Debug

无效值

为什么会出现 invalid value？在 Numpy 的计算中，只有两种可能会出现无效值：一是人为引入了无效值，这些无效值在参与计算的时候会报警；二是计算出现了无效值，例如除零等问题。我在优化过程中确保不会人为引入 np.inf，np.nan 这样的无效值，也确保了计算时不会出现除零问题，那这里又是为什么会出现无效值呢？很显然是开根号出现了负数，也就是 r_out_parallel 向量的长度的平方超过了 1。我们看一下这个函数的代码：

1def refract(self, normal: Vec3, etai_over_etat: float) -> Vec3:
2    cos_theta: float = -self @ normal
3    r_out_parallel: Vec3 = (self + normal*cos_theta) * etai_over_etat
4    r_out_prep: Vec3 = \
5        normal * (-np.sqrt(1 - r_out_parallel.length_squared()))
6    return r_out_parallel + r_out_prep

在这个函数中，我们会首先将单位长度的自己（入射光线）和单位长度的法线进行点乘并取负，获得 cosθ 的值（θ 为入射角）。然后我们根据折射定律，就可以求出折射光线的两个垂直分量的值。在求 r_out_parallel 的过程中，etai_over_etat 是一个根据折射率确定的定值，变量在于 self 和 normal，也就是入射光线和法线。如果两者都是确保为单位长度，那计算结果的长度绝对不会大于 1。

查看唯一调用了这个函数的 Dielectric 材质的代码，可以看到入射光线由这一行定义：

1unit_direction: Vec3 = r_in.direction().unit_vector()

非常明显，入射光线保证是单位长度。查看之前法线的计算代码：

1outward_normal: np.ndarray = (point - self.center.e) / self.radius

法线也确保是单位长度…？那这就非常奇怪了。

自然而然，我在 refract 函数加了输出 log 的代码，然后试图复现。由于太过随机，我只能开一个足够复杂的渲染，然后期待可能长达一个多小时的渲染中能够出现一次这个问题。结果证明，是法线可能出现长度大于 1 的情况。（而且不是误差的那种大一点点，是很明显的大很多；毕竟大一点点也不会触发报警）

法线长度

奇怪？明明我们的法线除以了球体的半径了，这应该会保证法线长度为单位长度的啊。我在这代码附近做了些测试，发现一般情况下，法线的长度可能会因为浮点误差稍稍大于 1，绝大部分的情况下长度都满足小于 1.003，这在意料之中，并不会造成问题。但是造成问题的就是在极稀有的情况下，这个长度会远超过 1，以下是记录到的几个异常向量的值（后面跟的是长度的平方）：

1[ 3.3827233  -1.2850341  -0.84744215] 13.812287
2[-0.08962154 -2.7663574   0.6211221 ] 8.046557
3[-1.1635733  -1.8047484   0.29541224] 4.698288
4[ 1.272459  -1.310974   1.5275431] 5.6711926

光看这个值看不出什么名堂，我一起记录了光线的入射点、方向、长度的信息：（第一条没有记录到）

Divergence 1

1normal: [-0.08962154 -2.7663574   0.6211221 ]
2length_squared: 8.046557
3ray.origin: [764.490478515625 -1161.707275390625 624.0230102539062]
4ray.direction: [-0.7604150772094727 1.1616960763931274 -0.6264601945877075]
5ray.t: 999.7055053710938

Divergence 2

1normal: [-1.1635733  -1.8047484   0.29541224]
2length_squared: 4.698288
3ray.origin: [-858.3599243164062 -920.5106201171875 -506.8523864746094]
4ray.direction: [0.8497467637062073 0.9204732179641724 0.5076704025268555]
5ray.t: 999.8658

Divergence 3

1normal: [ 1.272459  -1.310974   1.5275431]
2length_squared: 5.6711926
3ray.origin: [963.6956787109375 -784.3375854492188 157.41766357421875]
4ray.direction: [-0.9603748917579651 0.7843140959739685 -0.16341014206409454]
5ray.t: 999.9506225585938

首先能注意到的异状非常明显：光线的来源点非常可疑。这个场景在渲染的时候，地面是一个半径 1000 长度的大球，而这几个点非常明显都在球上，而且是很远很靠下的位置，渲染的时候光线应该不太能达到这些点。不过我们先假设能打到好了，计算本身应该不会有问题吧？在哪儿的点打过来都应该计算正确才对，但是这些光线算出来的法线就是会明显偏大。难道是误差？

np.float32

没错，我用了 np.float32。因为心想着单精度算起来更快，而且精度大概也够，双精度仿佛太 overkill，所以在设计之初选择了在系统里统一使用单精度。我们看一下计算光线与球体碰撞的代码：

1oc = r.origin() - self.center.e
2a = (r.direction()**2).sum(axis=1)
3half_b = (oc * r.direction()).sum(axis=1)
4c = (oc**2).sum(axis=1) - self.radius**2
5discriminant_list = half_b**2 - a*c

我们再看一看 Divergence 3 的计算中，这几个中间值：

1oc: [ 960.582   -784.5376   163.70758]
2a: 1.5641714
3half_b: -1564.5942
4c: 1565017.1
5discriminant: 0.25

震撼猫咪。根据 IEEE 754 标准，32 位浮点数用于表示有效数字的部分只有 23 bit，这个数字太大，以至于整数部分撑满了有效数字，小数部分完全不够表示了。于是这直接导致了最后计算 t 值时偏差过大，又导致计算得到的入射点偏差过大，最终导致法线的超大误差。

破案🎉！这种选择单精度浮点数的行为是懒惰的体现，是对万千先人对精度的追求的侮辱，这种人渣行为应当取缔。（但是真的是这样吗？）

“地球”

还有几个疑点我们没有解决：入射光线的源点，和光线的长度。为什么源点都在底下的那个大球（地球）上？为什么都在似乎完全照射不到的位置（我的模型完全不会考虑爱因斯坦老爷子讲的引力作用）？以及还有一个相当可疑的，为什么所有光线的长度都正好是 999.*？

我们先解决最容易的，999。整个光追器里，唯一一个值是 1000 的常量，就是地球的半径。这种奇妙的联系一定会怀疑吧？我尝试了一个很简单的关联测试，就是把地球的半径调整到 1200，然后再试着复现一次。

Divergence 4

1normal: [-2.146833   -4.778076    0.47927737]
2length_squared: 27.668612
3ray.origin: [-615.9678955078125 -2229.79443359375 10.366312980651855]
4ray.direction: [0.5196393132209778 1.858138084411621 -0.011417373083531857]
5ray.t: 1199.6087646484375

于是，这关联就板上钉钉实锤了。

到这份上，想要解释也不难了。想一下光线的源点和方向，就会发现这几个异常无一例外都是在地球表面很偏远的角落，并且法线指向球内部。我们看一看采用 lambertian 散射的磨砂材质的代码：

1class Lambertian(Material):
2    def __init__(self, a: Color) -> None:
3        self.albedo = a
4
5    def scatter(self, r_in: Ray, rec: HitRecord) -> Tuple[Ray, Color]:
6        scatter_direction: Vec3 = rec.normal + Vec3.random_unit_vector()
7        scattered = Ray(rec.p, scatter_direction)
8        attenuation = self.albedo
9        return scattered, attenuation

果然，我们没有考虑光线在球内部的情况。如果有光线进入了球内，这光线会继续按照漫反射的规则做随机发散。这也顺便解释了为什么异常光线长度都是地球半径：Lambertian 反射在这里的方向计算方式，是首先获得一个以法线末端为球心的单位球；然后在这球上按均匀分布随机取一点，则反射方向就是法线起始端（光线源点）指向这个球上的随机点。原书示意图如下，只不过这里的情况是单位球在大球内部。

这个单位球，显然和大球相似，相似比为 1：大球半径。因此最终的光线长度即为大球的半径。

Shadow Acne

光线长度的问题解决以后，最后一个问题就是：为什么会有光线进入地球内部？

在创建场景的时候，地球采用了磨砂材质，这是一种默认并不会允许光穿透过去的材质。从上一节可以看到，lambertian 散射的计算保证了反射光线一定在切平面之上，不可能进入球内。那这是从哪里来的光线？蓝色空间号从四维带进去的吗？

这里的突破口在主函数渲染函数 ray_color() 中的一行：

1rec_list: HitRecordList = world.hit(r, 0.001, np.inf)

原书 8.4 节介绍了这一方法，通过把光线最小距离设置为 0.001 以解决 Shadow Acne 的问题。这就导致了光线会突破限制进入地球内部。在计算光线和球体碰撞的时候，一元二次方程有解时往往会有两个解，我们使用了这样的逻辑去挑选最终的解：（代码为未并行化的版本）

1if discriminant > 0:
2    if t_min < t_0 < t_max:
3        t = t_0
4    elif t_min < t_1 < t_max:
5        t = t_1
6    else:
7        return None

当光线长度小于 0.001 时，它就会被抛弃。于是在地球和其他小球的接缝处，非常短的反射光线的第一个解被抛弃了，而第二个解就在地球的另一边。这也就是为什么，光线会偶然突破限制进入不应该进入的地球内部。

尾声

最终的修复也很简单，在 Lambertian 材质的计算中加入一个光线位置的判断即可。我们本来就在光线的击中记录中加入了在球内还是球外的记录项，因此非常方便。

最后附上一张 debug 时渲染的测试图作为尾声：