[Maxim07]中光线与三角形求交算法的推导

"Ray-Triangle Intersection Algorithm for Modern CPU Architectures" Maxim Shevtsov, Alexei Soupikov and Alexander Kapustin (2007) 中的公式推导时的约定与我通常使用的约定不一致，并且叙述并不严格，故自己推导如下：

约定

以下公式均定义在右手系中。

* 表示数乘或点乘（内积），x表示叉乘（外积）

光线（射线）的参数表示为 x = o + d * t (t > 0) (1)

其中x为光线上任意一点，o为光线起点，d为光线方向，t为参数。

三角形的参数表示为 x = (1 - u - v) * p + u * p1 + v * p2 (0 <= u <= 1, 0 <= v <= 1) (2)

x为三角形内的任一点，p, p1, p2为三角形的三个顶点位置，u, v为参数，（1 - u - v, u, v)称为三角形的重心坐标（barycentric coordinates）。

为了方便，令e1 = p1 - p, e2 = p2 - p

我们约定p, p1, p2按顺序构成逆时针的方向为正面，即n = e1 x e2

 

推导

光线与三角形的交点可以表示为线性方程 (1 - u - v) * p + u * p1 + v * p2 = o + d * t (3)

即 e1 * u + e2 * v - d * t = o - p (4)

设K = -d, L = o - p, 则原式变为 e1 * u + e2 * v + K * t = L (5)

利用Cramer's rule解该线性方程得

t = dett / det, u = detu / det, v = detv / det (6)

其中

det = K * (e1 x e2) = - d * n (7)

dett = L * (e1 x e2) = (o - p) * n (8)

detu = K * (L x e2) = -d * [(o - p) x e2] (9)

detv = K * (e1 x L) = -d * [e1 x (o - p)] (10)

其实到这里，已经求出了判断交点所需的信息，当然基于以上公式，我们还会预计算一些信息，储存在三角形中，减少求交所需的计算量。

事实上，我们可以在公式(6)的分子分母同除以一个数，所得的值不变，利用这一点，我们同除以nw，w是n的最大分量的轴向，令u, v为另外两个轴向，且u < v。

故 det / nw = -(du * nu + dv * nv + dw) (11)

dett / nw = (ou * nu + ov * nv + ow) - (pu * nu + pv * nv + pw) (12)

detu / nw = d * [(p - o) x e2] = d * [p x e2 - o x e2] = du * (pv * e2w - pw * e2v - ov * e2w + ow * e2v) + dv * (pw * e2u - pu * e2w - ow * e2u + ou * e2w) + dw * (pu * e2v - pv * e2u - ou * e2v + ov * e2u) (13)

利用e1 * n = 0, e2 * n = 0我们可以得到

e1w = - (e1u * nu + e1v * nv) (14)

e2w = - (e2u * nu + e2v * nv)

利用公式(11)，我们可以得到

dw = -(du * nu + dv * nv + det / nw) (15)

将公式(14), (15)代入(13)中消去e2w, dw并整理出e2u, e2v得

detu / nw = e2v * [du * (ou - pu) * nu + du * (ov - pv) * nv + dw * (ow - pw) - (pu - ou) * det / nw] - e2u * [dv * (ou - pu) * nu + dv * (ov - pv) * nv + dv * (ow - pw) - (pv - ov) * det / nw] (16)

令Du = du * dett / nw - (pu - ou) * det, Dv = dv * dett / nw - (pv - ov) * det代入(16)得

detu / nw = e2v * Du - e2u * Dv (17)

同理，实际上我们不需要计算，只需要观察公式(9)与公式(10)的区别，把e2换成e1并加一负号得

detv / nw = e1u * Dv - e1v * Du (18)

至此，我们已经推导出了整个算法中需要用到的公式，下节给出算法的代码。

实现

实现中我们为每个三角形预计算并储存如下数据结构，正好是48 Byte大小，可以对齐到16 Byte边界：

struct TriAccel {
    uint w;
    uint nw;
    float np;        
    float nu;
    float nv;
    float pu;
    float pv;
    float e1u;
    float e2v;
    float e2u;
    float e1v;
    uint tri_id;
};

预计算三角形代码如下，为了处理背面裁剪增加了一些代码，不多介绍了：

TriAccel::TriAccel(const Vector3f & v1_pos, const Vector3f & v2_pos, const Vector3f & v3_pos)
{
    Vector3f e1;
    Vector3f e2;
    Vector3f n;
    Vector3f absn;
    float maxVal;
    float inv_nw;
    uint u, v;

    sub(e1, v2_pos, v1_pos);
    sub(e2, v3_pos, v1_pos);
    cross(n, e1, e2);
    absn.x = absf(n.x);
    absn.y = absf(n.y);
    absn.z = absf(n.z);

    if (absn.x > absn.y) {
        u = Y_AXIS;
        v = Z_AXIS;
        w = X_AXIS;
        maxVal = absn.x;
    } else {
        u = Z_AXIS;
        v = X_AXIS;
        w = Y_AXIS;
        maxVal = absn.y;
    }
    if (absn.z > maxVal) {
        u = X_AXIS;
        v = Y_AXIS;
        w = Z_AXIS;
    }

    inv_nw = 1.0f / n.arr[w];
    mulvf(n, inv_nw);

    nu = n.arr[u];
    nv = n.arr[v];

    pu = v1_pos.arr[u];
    pv = v1_pos.arr[v];

    np = n.arr[u] * v1_pos.arr[u] + n.arr[v] * v1_pos.arr[v] + v1_pos.arr[w];

    e1u = (e1.arr[u] * inv_nw);
    e1v = (e1.arr[v] * inv_nw);

    e2u = (e2.arr[u] * inv_nw);
    e2v = (e2.arr[v] * inv_nw);

    if (inv_nw >= 0.0f)
    {
        nw = 0x00000000;
    }
    else
    {
        nw = 0x80000000;
    }
}

最后是光线与三角形求交代码：

 

#define FLOAT_TO_DWORD(x) (*(unsigned long*)(&(x)))

#define FLOAT_BITWISE_XOR(dest, lhs, rhs) \
    FLOAT_TO_DWORD(dest) = (FLOAT_TO_DWORD(lhs)) ^ (FLOAT_TO_DWORD(rhs))

#define FLOAT_BITWISE_OR(dest, lhs, rhs) \
    FLOAT_TO_DWORD(dest) = (FLOAT_TO_DWORD(lhs)) | (FLOAT_TO_DWORD(rhs))

#define FLOAT_MUL_NEGATIVE(lhs, rhs) \
    ((FLOAT_TO_DWORD(lhs) ^ FLOAT_TO_DWORD(rhs)) & 0x80000000) == 0x80000000
 
static const uint sAxisTable[3][2] = {
    {Y_AXIS, Z_AXIS}, 
    {Z_AXIS, X_AXIS}, 
    {X_AXIS, Y_AXIS}, 
};

bool TriAccel::intersect(const Vector3f & ray_src, const Vector3f & ray_dir, const float t_near, const float t_far, float & t, Vector3f & bary)
{
    float det;
    float dett;
    float Du;
    float Dv;
    float detu;
    float detv;
    float tempdet0;
    float tempdet1;
    uint u, v;

    u = sAxisTable[w][0];
    v = sAxisTable[w][1];

    det = -(ray_dir.arr[u] * nu + ray_dir.arr[v] * nv + ray_dir.arr[w]);

    if (FLOAT_MUL_NEGATIVE(det, nw))
    {
        return false;
    }

    dett = (ray_src.arr[u] * nu + ray_src.arr[v] * nv + ray_src.arr[w]) - np;

    if (FLOAT_MUL_NEGATIVE(det, dett))
    {
        return false;
    }

    Du = ray_dir.arr[u] * dett - (pu - ray_src.arr[u]) * det;
    Dv = ray_dir.arr[v] * dett - (pv - ray_src.arr[v]) * det;
    detu = e2v * Du - e2u * Dv;
    detv = e1u * Dv - e1v * Du;

    tempdet0 = det - detu - detv;

    FLOAT_BITWISE_XOR(tempdet0, tempdet0, detu);
    FLOAT_BITWISE_XOR(tempdet1, detv, detu);
    FLOAT_BITWISE_OR(tempdet0, tempdet0, tempdet1);

    if ((FLOAT_TO_DWORD(tempdet0) & 0x80000000) == 0x80000000)
    {
        return false;
    }

    float inv_det = 1.0f / det;

    t = dett * inv_det;
    if (t > t_far || t < t_near)
    {
        return false;
    }

    bary.y = detu * inv_det;
    bary.z = detv * inv_det;
    bary.x = 1.0f - (bary.y + bary.z);

    return true;
}

到此该睡觉去鸟……文中如有错漏之处，欢迎指正，谢谢！