在光线追踪中避免自相交的方法
这是我阅读 Ray Tracing Gem 的一篇笔记，《避免自相交的快速可靠的方法》是 Ray Tracing Gem 的第六章。这篇博客文章主要是为了记录和解释原文末尾给出的魔法一般的代码。
浮点误差
数学上给出了很多不同的求解光线与几何体相交的方法，它们都是等价的，但对于计算机而言并不一定。由于浮点数带来的浮点误差，不同方法带来的结果可能天差地别。
通常来讲，浮点数越大，其相对应的误差也就越大，这不难理解。IEEE754 单精度浮点数由 1 位符号位、8 位阶码以及 23 位尾数组成。浮点数的数量级主要是由阶码决定的，因为阶码的增减会以指数级别影响浮点数的大小。所以，如果一个浮点数的数量级越大，其阶码就越大，那么尾数的舍入的影响也就越大。
避免自相交的方法
图元 ID 排除法
可以使用ID显式地标记已经相交过的图元。
优点：
- 不需要任何参数；
- 满足缩放不变性；
- 不会略过临近的几何体。
缺点：
- 若交点落在两个图元的交线处，或者出射光纤与平面夹角很小，则仍然会导致自相交；
- 无法处理重复或者重叠的几何体；
- 仅可用于处理平面图元，因为非平面图元可能产生有效自相交。
限制光线区间
可以限制光线的有效相交区间。这样，当光线距离相交平面过近时，相交会被排除。
优点：
- 几乎没有额外开销
缺点：
- 非常不可靠，对于不同的场景会有不同的结果；
- 可能导致小夹角下的自相交；
- 可能导致略过相邻平面；
沿几何法向量进行自适应偏移
这是本文要介绍的重点，来自NVIDIA的魔法。
前面在浮点误差处提到过，交点的误差与光源到交点的距离有关，所以自适应偏移的距离也应当与这个距离有关。使用固定偏移量\(\epsilon\)显然无法自动适应距离的变化。因此，他们统计了光源到交点的距离与误差的关系：
容易发现，当相交点与光源的距离足够大时，距离和误差的数量级之间基本上呈现线性关系。当相交点与光源的距离过小时，则使用固定的距离进行偏移。Ray Tracing Gem给出的代码实现如下：
constexpr float origin() { return 1.0f / 32.0f; } // 距离界限
constexpr float float_scale() { return 1.0f / 65536.0f; } // 用于固定偏移
constexpr float int_scale() { return 256.0f; } // 数量级差距
float3 offset_ray(const float3 p, const float3 n)
{
int3 of_i(int_scale() * n.x, int_scale() * n.y, int_scale() * n.z);
float3 p_i(
int_as_float(float_as_int(p.x) + ((p.x < 0) ? -of_i.x : of_i.x)),
int_as_float(float_as_int(p.y) + ((p.y < 0) ? -of_i.y : of_i.y)),
int_as_float(float_as_int(p.z) + ((p.z < 0) ? -of_i.z : of_i.z)));
return float3(fabsf(p.x) < origin() ? p.x + float_scale() * n.x : p_i.x,
fabsf(p.y) < origin() ? p.y + float_scale() * n.y : p_i.y,
fabsf(p.z) < origin() ? p.z + float_scale() * n.z : p_i.z);
}
这里面用了魔法般的int_as_float和float_as_int，导致这段代码非常令人迷惑。int_as_float用于按照float的格式解析这段int数据，相当于*reinterpret_cast<float *>(&int_value)，float_as_int同理。
为了理解上面的代码，我们对其进行一定的简化——我们将向量换成float，暂且去掉最后一部分，假设n是normalized vector，且为正数：
static constexpr float origin = 1.0f / 32.0f;
static constexpr float float_scale = 1.0f / 65536.f;
static constexpr float int_scale = 256.0f;
static int float_as_int(float value) {
union {
float f;
int i;
} v;
v.f = value;
return v.i;
}
static float int_as_float(int value) {
union {
float f;
int i;
} v;
v.i = value;
return v.f;
}
float offset_ray(const float point, const float normal) {
int off_i = static_cast<int>(int_scale * normal);
float p_i = int_as_float(float_as_int(point) + off_i);
return p_i;
}
首先，normal是正则化后的向量，因此off_i一定不会超过256，也就是说off_i一定不超过0x0000 0100。
然后，我们考虑第二行的float_as_int(point) + off_i。在IEEE754单精度浮点数中，低23位表示尾数。因此，float_as_int(point) + off_i等价于point的尾数加off_i（因为off_i只有低12位有效，其他部分全是0）。而浮点数的表示为\((-1)^S \times 2^{T - 127} \times M\)，其中\(S\)表示符号位，\(T\)表示阶码，\(M\)表示尾数。因此，直接对尾数进行加减法会自动放大至\(2^{T-127}\)倍，也就是说会按照原来浮点数的数量级进行缩放。而上面的代码就是利用浮点数的这一特性来正确处理偏移量与距离的线性关系的。
至于最后一段代码：
return float3(fabsf(p.x) < origin() ? p.x + float_scale() * n.x : p_i.x,
fabsf(p.y) < origin() ? p.y + float_scale() * n.y : p_i.y,
fabsf(p.z) < origin() ? p.z + float_scale() * n.z : p_i.z);
则是特殊判断交点到光源的距离是否过近，如果过近则直接应用常数偏移，没有什么特殊的地方，不再详述。
存在的问题
这个方法仍然不能完美处理交点沿着几何法向量偏移后跳过三角面的情况，我们仍然能够构造特殊的情况使之失效。不过，在实际应用中，这个方法出现失效的情况相比前面几种方法要小得多。