射线—方盒相交计算 Ray-Box Intersection

最新推荐文章于 2023-03-26 21:28:55 发布

Mr_Darkside

最新推荐文章于 2023-03-26 21:28:55 发布

阅读量4.2k

点赞数 1

分类专栏： 3D数学文章标签：数学 AABB相交测试 C++

3D数学专栏收录该内容

0 篇文章 0 订阅

订阅专栏

说在前面：

我的第一篇博客文章，先翻译一篇外网的文章，顺便锻炼一下翻译能力。

英语水平有限，有能力的读者建议看原文。

译文：

射线-方盒相交计算

下面的例子我们将假设这个方盒对齐于我们的坐标系统的坐标轴。这样的方盒被称为“坐标轴对齐方盒”，因为这样的方盒通常作为一个边界包围盒（或者边界包围体），所以也更通常的称为AABB盒(axis-aligned bounding box，坐标轴对齐边界包围盒 / 轴对齐有界盒)。其他的那些不与笛卡尔坐标系的坐标轴指向一致的方盒被称为OBB盒(oriented bounding box，方向性边界包围盒 / 有向有界盒)。

计算一条射线与AABB的交点是非常简单的。我们唯一所要记住的一点是，一条直线可以被下面这样的解析式所定义：

y = mx + b

从数学上讲，m 项叫做坡度（或者梯度），它负责表达直线的指向。b 项对应于直线与y轴的交点。

这条直线也可以用下面这个方程表示：

O + Rt

O对应射线的起点，R是它的方向，参数t可以是任何实数(负数，正数或者零)。只要随意调节方程中的t，我们就可以创造被这条用位置与方向定义的射线上的任何一个点。这个射线方程与一般直线方程十分相似。实际上当你用b和m替换掉O和R时它们是一样的。为了表示一个坐标轴对齐的边界包围盒，我们所需的仅是两个用来表示方盒的最大和最小范围的点(在代码中被称为bounds)。

template<typename T>
class Box3
{
public:
    Box3(Vec3<T> vmin, Vec3<T> vmax)
    {
        bounds[0] = vmin;
        bounds[1] = vmax;
    }
    Vec3<T> bounds[2];
};

这个立体的bounds定义了一系列平行于坐标轴的直线，我们也可以使用直线方程来表示它们。举个例子，边界包围盒的最小范围关于x分量的直线方程可以被这样写：

y = B0x

为了找到射线与这条直线的交点可以这么写：

Ox + tRx = B0x (方程1)

这可以移项：

t0x = (B0x - Ox) / Rx (方程2)

这种计算边界包围盒的最小范围的x分量的方法也可以相似的用于计算t1x。注意当t的值是负的时候，方盒在射线的后面。

最终如果我们在这个过程的结尾对y和z分量应用相同的技巧，将得到一组六个值，指明射线与方盒平行于x，y和z轴的平面的相交位置。

t0x = (B0x - Ox) / Rx

t1x = (B1x - Ox) / Rx

t0y = (B0y - Oy) / Ry

t1y = (B1y - Oy) / Ry

t0z = (B0z - Oz) / Rz

t1z = (B1z - Oz) / Rz

注意如果射线平行于坐标轴，那么它将不会与在这个坐标轴上的边界包围盒的面相交(在这种情况下，这条线的方程会缩减到常数b并且方程1无解)。这种情形下的问题是要找出六个值中的哪一个相对应于射线与方盒的交点(是否射线与方盒相交了。目前为止，我们只是找到了射线与被一个立方体的每一个面所定义的平面的交点)。可以对每一个计算出来的t进行一个简单的测试来简单地找到答案。从插图2中可以看出这背后的逻辑是非常明显的(这个例子中我们只考虑2D下的情况)。这条射线先与那个被方盒的最小范围所定义的平面在两个位置相交：t0x和t0y。然而，与这些平面的相交并不一定意味着这些交点在立方体上(如果它们不在立方体上，很明显射线不与方盒相交)。从数学角度，我们可以通过比较它们的值来找出其中哪一个位于立方体上：我们只需要简单地找出那个t值更大一些的交点。用伪代码我们这样写：

t0x = (B0x - Ox) / Rx
t0y = (B0y - Oy) / Ry
tmin = (t0x > t0y) ? t0x : t0y

找到射线与方盒第二个交点的过程是类似的，不过在那种情况下，我们将使用被方盒最大范围所定义的平面计算出来的t值，并且选择t值最小的那个点。

t1x = (B1x - Ox) / Rx
t1y = (B1y - Oy) / Ry
tmax = (t1x < t1y) ? t1x : t1y

然而，射线不总是与方盒相交。我们在插图3中展示了几种射线错过方盒的情况。这些情况也都可以用简单比较t的值的方法来分清。正如你在图中所能看到的，当t0x大于t1y并且t0y大于t1x的时候，射线错过了方盒。

让我们把这个测试加到我们的代码中：

if (t0x > t1y || t0y > t1x)
   return false

最终我们可以把这种技巧拓展到3D情形，通过为z分量计算t值，并与目前我们已为x和y分量计算的值所比较：

t0z = (B0z - Oz) / Rz
t1z = (B1z - Oz) / Rz
if (tmin > t1z || t0z > tmax) return false
if (t0z > tmin) tmin = t0z
if (t1z < tmax) tmax = t1z

到这一步，我们已经确定了射线会和方盒相交。然而，你可能需要为射线设置一个tmin和tmax值，这取决于你所使用的实现方式。在返回(return)射线与方盒的交点之前，你可能需要使用事先为射线设定的tmin和tmax值来检查你在射线-方盒交点处算出来的最大t值和最小t值：

if (tmin > ray.tmax || tmax < ray.tmin)
   return false

如果我们通过了测试，那么就可以使用来自于函数的tmin和tmax来设置最小和最大t值了。这些都对应于射线与方盒交点处t的参数值：

if (tmin > ray.tmin) ray.tmin = tmin
if (tmax < ray.tmax) ray.tmax = tmax

这里是C++中完整的实现方法(min和max是边界包围盒的最小和最大范围)：

bool intersect(const Ray<T> &r)
{
    T tmin = (min.x - r.orig.x) / r.dir.x;
    T tmax = (max.x - r.orig.x) / r.dir.x;
    if (tmin > tmax) swap(tmin, tmax);
    T tymin = (min.y - r.orig.y) / r.dir.y;
    T tymax = (max.y - r.orig.y) / r.dir.y;
    if (tymin > tymax) swap(tymin, tymax);
    if ((tmin > tymax) || (tymin > tmax))
        return false;
    if (tymin > tmin)
        tmin = tymin;
    if (tymax < tmax)
        tmax = tymax;
    T tzmin = (min.z - r.orig.z) / r.dir.z;
    T tzmax = (max.z - r.orig.z) / r.dir.z;
    if (tzmin > tzmax) swap(tzmin, tzmax);
    if ((tmin > tzmax) || (tzmin > tmax))
        return false;
    if (tzmin > tmin)
        tmin = tzmin;
    if (tzmax < tmax)
        tmax = tzmax;
    if ((tmin > r.tmax) || (tmax < r.tmin)) return false;
    if (r.tmin < tmin) r.tmin = tmin;
    if (r.tmax > tmax) r.tmax = tmax;
    return true;
}

注意tmin可能大于tmax，这取决于射线的方向。考虑到之后我们所要进行的测试的所有逻辑，都依赖于t0始终比t1小，所以当t1比t0小的时候我们需要交换它们的值(5,8和17行)。

优化代码

我们可以对前面的代码做很多改进，使它不仅变得更快而且更健壮。首先我们可以用下面的代码替换交换操作：

if (ray.dir.x >= 0) {
    tmin = (min.x - r.orig.x) / ray.dir.x;
    tmax = (max.x - r.orig.x) / ray.dir.x;
}
else {
    tmin = (max.x - r.orig.x) / ray.dir.x;
    tmax = (min.x - r.orig.x) / ray.dir.x;
}

然而这个代码有一个问题。当射线的方向是0的时候会造成被0除。这应该由编译器返回+∞来恰当解决。当射线的方向是-0的时候问题发生了。当射线的方向是负的时候，我们希望if语句的第二块代码被执行，但是万一当射线方向是-0时，反而第一块代码会被执行，因为在IEEE浮点格式里，测试-0 = 0会返回真。值tmin/tmax会被设置成+∞和-∞(而不是-∞/+∞),并且代码不会成功检测到相交。我们可以用射线方向的倒数替换射线方向来简单地解决这个问题:

Vec3<T> invdir = 1 / ray.dir;
if (invdir.x >= 0) {
    tmin = (min.x - r.orig.x) * invdir.x;
    tmax = (max.x - r.orig.x) * invdir.x;
}
else {
    tmin = (max.x - r.orig.x) * invdir.x;
    tmax = (min.x - r.orig.x) * invdir.x;
}

当射线的方向是-0的时候方向的倒数就应该是-∞，并且第二行的测试将如我们所愿地返回假。为了应对射线与许多方盒检测相交的情况，我们可以在射线的构造函数中预计算射线方向的倒数，然后在之后的相交中重复使用来节省一些时间。我们也可以预计算射线方向的符号，从而以一种更加紧凑的形式来写这个方法。这里是射线-方盒相交方法的最终版本：

template<typename T>
class Ray
{
public:
    Ray(Vec3<T> orig, Vec3<T> dir) : orig(orig), dir(dir), tmin(T(0)), tmax(std::numeric_limits<T>::max())
    {
        invdir = T(1) / dir;
        sign[0] = (invdir.x < 0);
        sign[1] = (invdir.y < 0);
        sign[2] = (invdir.z < 0);
    }
    Vec3<T> orig, dir;       /// ray orig and dir
    mutable T tmin, tmax;    /// ray min and max distances
    Vec3<T> invdir;
    int sign[3];
};
 
bool intersect(const Ray<T> &r) const
{
    T tmin, tmax, tymin, tymax, tzmin, tzmax;
    tmin = (bounds[r.sign[0]].x - r.orig.x) * r.invdir.x;
    tmax = (bounds[1-r.sign[0]].x - r.orig.x) * r.invdir.x;
    tymin = (bounds[r.sign[1]].y - r.orig.y) * r.invdir.y;
    tymax = (bounds[1-r.sign[1]].y - r.orig.y) * r.invdir.y;
    if ((tmin > tymax) || (tymin > tmax))
        return false;
    if (tymin > tmin)
        tmin = tymin;
    if (tymax < tmax)
        tmax = tymax;
    tzmin = (bounds[r.sign[2]].z - r.orig.z) * r.invdir.z;
    tzmax = (bounds[1-r.sign[2]].z - r.orig.z) * r.invdir.z;
    if ((tmin > tzmax) || (tzmin > tmax))
        return false;
    if (tzmin > tmin)
        tmin = tzmin;
    if (tzmax < tmax)
        tmax = tzmax;
    if (tmin > r.tmin) r.tmin = tmin;
    if (tmax < r.tmax) r.tmax = tmax;
    return true;
}

这个优化过的版本在我们的机器上运行比第一个版本快1.5倍(速度的增加可能依赖于编译器)。

这些代码返回射线与方盒的交点，它可能在射线原点的前面或者后面。举个例子，如果射线的原点在方盒内部(就像左边的图片)，将会有两个交点：一个在射线的前面，一个在后面。我们知道当t的值是负的时候，交点在射线原点的后面。当t是正的时候，交点在射线原点的前面。如果你的算法对于t小于0的交点不感兴趣，那么你将不得不在返回射线-方盒交点时，小心地处理这些情况。(这常常是bug的来源)