11 - 射线和三角形的相交性检测

1. 步骤

• 计算射线和包含该三角形的平面的交点；
• 通过计算交点的重心坐标，判断其是否在三角形中。

2. 优化

• 在检测中尽可能早地返回负值，这称作 “提前结束（early out）” ；
• 尽可能地延迟 “昂贵” 的数学运算，如除法。
  • 若并不需要 “昂贵” 运算的结果，比如说遇到了 “提前结束” 的情况，那么执行这些运算的时间就白白浪费了。
  • 它给了编译器更多的空间以利用现代处理器的指令管道的优点。若指令（如除法）有很长的延迟，编译器就能向前查看并预先生成开始除法运算的代码。在准备除法运算时，它产生执行其他测试的代码（可能导致提前结束）。所以，在执行期间，如果确实需要除法运算的结果，该结果可能已经被计算出来，或至少已经部分被计算出来了。
• 只检测与三角形正面的相交，这几乎可以节省一半的检测时间。

3. 代码

// 射线与三角形的相交性检测
float rayTriangleIntersect(
    const Vector3 &rayOrg,      // 射线起点
    const Vector3 &rayDelta,    // 射线长度和方向
    const Vector3 &p0,          // 三角形顶点
    const Vector3 &p1,          // .
    const Vector3 &p2,          // .
    float minT                  // 目前为止最近的点，从 1.0 开始
){
    // 如果没有相交就返回这个大数
    const float kNoIntersection = 1e30f;
    // 计算顺时针的边向量
    Vector3 e1 = p1 - p0;
    Vector3 e2 = p2 - p1;
    // 计算表面法向量（不需要正则化）
    Vector3 n = crossProduct(e1, e2);
    // 计算倾斜角，表示靠近三角形“正面”的程度
    float dot = n * rayDelta;
    // 检查射线平行于三角形或未指向三角形正面的情况
    // 注意这将拒绝退化三角形和射线，这里的编码方式非常特殊，NAN 将不能通过检测（当涉及到 NAN 时，与 dot >= 0.0f 时的表现不同）
    if(!(dot < 0.0f)){
        return kNoIntersection;
    }
    // 计算平面方程的 d 值，在右边使用带 d 的平面方程
    // 
    // Ax + By + Cz = d
    float d = n * p0;
    // 计算和包含三角形的平面的参数交点
    float t = d - n * rayOrg;
    // 检查射线起点是否在多边形的背面
    // 再次，我们显示这个检查以便确保 NANs
    if(!(t <= 0.0f)){
        return kNoIntersection;
    }
    // 交点已经找到（或射线并未到达平面）？
    // 
    // 因为 dot < 0;
    // 
    // t/dot > minT
    // 与
    // t < dot * minT
    // 是相同的
    // （之后我们将其倒置以便进行 NAN 检查。。。）
    if(!(t >= dot * minT)){
        return kNoIntersection;
    }
    // 射线和平面相交，计算实际的交点
    t /= dot;
    assert(t >= 0.0f);
    assert(t <= minT);
    // 计算 3D 交点
    Vector3 p = rayOrg + rayDelta * t;
    // 找到最主要的轴，选择投影平面
    float u0, u1, u2;
    float v0, v1, v2;
    if(fabs(n.x) > fabs(n.y)){
        // 以 x 轴为主轴，对平面 yoz 进行投影
        if(fabs(n.x) > fabs(n.z)){
            u0 = p.y - p0.y;
            u1 = p1.y - p0.y;
            u2 = p2.y - p0.y;
            v0 = p.z - p0.z;
            v1 = p1.z - p0.z;
            v2 = p2.z - p0.z;
        // 以 z 轴为主轴，对平面 xoy 进行投影
        }else{
            u0 = p.x - p0.x;
            u1 = p1.x - p0.x;
            u2 = p2.x - p0.x;
            v0 = p.y - p0.y;
            v1 = p1.y - p0.y;
            v2 = p2.y - p0.y;
        }
    }else{
        // 以 x 轴为主轴，对平面 yoz 进行投影
        if(fabs(n.y) > fabs(n.z)){
            u0 = p.x - p0.x;
            u1 = p1.x - p0.x;
            u2 = p2.x - p0.x;
            v0 = p.z - p0.z;
            v1 = p1.z - p0.z;
            v2 = p2.z - p0.z;
        // 以 z 轴为主轴，对平面 xoy 进行投影
        }else{
            u0 = p.x - p0.x;
            u1 = p1.x - p0.x;
            u2 = p2.x - p0.x;
            v0 = p.y - p0.y;
            v1 = p1.y - p0.y;
            v2 = p2.y - p0.y;
        }
    }
    // 计算分母，检查其有效性
    float temp = u1 * v2 - v1 * u2;
    if(!(temp != 0.0f)){
        return kNoIntersection;
    }
    // 计算重心坐标，每一步都检查边界条件
    float alpha = (u0 * v2 - v0 * u2) * temp;
    if(!(alpha >= 0.0f)){
        return kNoIntersection;
    }
    float beta = (u1 * v0 - v1 * u0) * temp;
    if(!(beta >= 0.0f)){
        return kNoIntersection;
    }
    float gamma = 1.0f - alpha - beta;
    if(!(gamma >= 0.0f)){
        return kNoIntersection;
    }
    // 返回参数交点
    return t;
}