射线与三角形相交算法(源于DirectX SDK)

算法来自：DirectX SDK, Picking程序示例, Pick.cpp文件。
此算法速度快、精度也高。

算法原理¶

射线参数方程¶

射线的参数方程如下，其中，R为射线终点，O是射线的起点，D是射线的方向：

\[ R = O + Dt \]

三角形参数方程¶

三角形的参数方程如下，其中

\(P\)表示三角形上的任意一点
\(V_0\)，\(V_1\)和\(V_2\)是三角形的三个点
\(u\)，\(v\)是\(V_1\)和\(V_2\)的权重，\(1-u-v\)是\(V_0\)的权重
满足\(u>=0,v>=0,v+v<=1\)

\(P = (1-u-v)V_0+uV_1+vV_2\)

确切的说，上面的方程是三角形及其内部所有点的方程

因为三角形内任意一点都可以理解为从顶点\(V_0\)开始
沿着边\(V_1V_2\)移动一段距离
然后再沿着边\(V_0V_2\)移动一段距离
然后求他们的和向量
至于移动多大距离，就是由参数\(u\)和\(v\)控制的

联立方程组¶

于是，求射线与三角形的交点也就变成了解下面这个方程

其中\(t\)，\(u\)，\(v\)是未知数，其他都是已知的

\(O+Dt=(1-u-v)V_0+uV_1+vV_2\)

移项并整理，将t,u,v提取出来作为未知数，得到下面的线性方程组

\[\begin{bmatrix} -D&V_1-V_0&V_2-V_0 \end{bmatrix} \begin{bmatrix} t\\u\\v \end{bmatrix} = O-V_0 \]

解方程组¶

现在开始解这个方程组，这里要用到两个知识点，一是克莱姆法则，二是向量的混合积。

令E1 = V1 - V0，E2 = V2 - V0，T = O - V0上式可以改写成

\[ \begin{bmatrix} -D&E_1&E_2 \end{bmatrix} \begin{bmatrix} t\\u\\v \end{bmatrix} = T \]

克莱姆法则¶

根据克莱姆法则，可得到t,u,v的解分别是

\[t= \cfrac{1}{\begin{vmatrix}-D&E_1&E_2\end{vmatrix}} \begin{vmatrix} T&E_1&E_2 \end{vmatrix} \]

\[u= \cfrac{1}{\begin{vmatrix} -D&E_1&E_2 \end{vmatrix}} \begin{vmatrix} -D&T&E_2 \end{vmatrix} \]

\[v= \cfrac{1}{\begin{vmatrix} -D&E_1&E_2 \end{vmatrix}} \begin{vmatrix} -D&E_1&T \end{vmatrix} \]

将这三个解联合起来写就是

\[\begin{bmatrix} t\\u\\v \end{bmatrix} = \cfrac{1}{\begin{vmatrix} -D&E_1&E_2 \end{vmatrix}} \begin{vmatrix} T &E_1 &E_2 \\ -D &T &E_2 \\ -D &E_1 &T \\ \end{vmatrix} \]

向量混合积¶

根据混合积公式

\[\begin{vmatrix} a&b&c \end{vmatrix} =a \times b \bullet c =-a \times c \bullet b \]

上式可以改写成

\[\begin{bmatrix} t\\u\\v \end{bmatrix} = \cfrac{1}{\begin{bmatrix} D \times E_2 \bullet E_1 \end{bmatrix}} \begin{bmatrix} T \times E_1 \bullet E_2 \\ D \times E_2 \bullet T \\ T \times E_1 \bullet D \\ \end{bmatrix} \]

令

\[ P=D \times E_2 \\ Q=T \times E_1 \]

最终的计算公式¶

得到最终的公式，这便是下面代码中用到的最终公式了，之所以提炼出P和Q是为了避免重复计算

\[\begin{bmatrix} t\\u\\v \end{bmatrix} = \cfrac{1}{\begin{vmatrix} P \bullet E_1 \end{vmatrix}} \begin{vmatrix} Q \bullet E_2 \\ P \bullet T \\ Q \bullet D \\ \end{vmatrix} \]

代码¶

bool intersects(const Ray& ray, const Triangle& triangle, 
                Vector3D* intersectionPnt = nullptr)
{
    const auto& O = ray.getStart();       
    const auto& dir = ray.getDirection(); 
    const auto& A = triangle.A();         
    const auto& B = triangle.B();         
    const auto& C = triangle.C();         

    // E1
    auto E1_AB = B - A;

    // E2
    auto E2_AC = C - A;

    // P
    auto P = dir ^ E2_AC;

    // determinant
    auto det = E1_AB * P;

    // keep det > 0, modify T accordingly
    Vector3D T;
    if (det > 0)
    {
        T = O - A;
    }
    else
    {
        T = A - O;
        det = -det;
    }

    // If determinant is near zero, ray lies in plane of triangle
    //射线与三角形所在平面平行
    if (det < 0.0001f)
    {
        return false;
    }

    // Calculate u and make sure u <= 1
    auto u = T * P; // barycentric coordinate of intersection
    if (u < 0.0 || u > det)
    {
        return false;
    }

    // Q
    auto Q = T ^ E1_AB;

    // Calculate v and make sure u + v <= 1
    auto v = dir * Q; // barycentric coordinate of intersection
    if (v < 0.0 || u + v > det)
    {
        return false;
    }

    // Calculate t, scale parameters, ray intersects triangle
    auto t = E2_AC * Q; // weight of the intersection for the ray

    double inv_det = 1.0f / det;
    t *= inv_det;
    u *= inv_det;
    v *= inv_det;

    if (t < 0)
    {
        //射线反方向上的交点
        return false;
    }

    if (intersectionPnt)
    {
        *intersectionPnt = O + dir * t;
    }
    return true;
}