"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边界：
三角形数据结构
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;
};
预计算三角形代码如下，为了处理背面裁剪增加了一些代码，不多介绍了：
预计算三角形
{
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_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;
}
到此该睡觉去鸟……文中如有错漏之处，欢迎指正，谢谢！