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

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

算法原理¶

射线参数方程¶

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

R = O + D t

三角形参数方程¶

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

$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} + u V_{1} + v V_{2}$

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

因为三角形内任意一点都可以理解为从顶点 $V_{0}$ 开始
沿着边 $V_{1} V_{2}$ 移动一段距离
然后再沿着边 $V_{0} V_{2}$ 移动一段距离
然后求他们的和向量
至于移动多大距离，就是由参数 $u$ 和 $v$ 控制的

联立方程组¶

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

其中 $t$ ， $u$ ， $v$ 是未知数，其他都是已知的

$O + D t = (1 - u - v) V_{0} + u V_{1} + v V_{2}$

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

[\begin{matrix} - D & V_{1} - V_{0} & V_{2} - V_{0} \end{matrix}] [\begin{matrix} t \\ u \\ v \end{matrix}] = O - V_{0}

解方程组¶

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

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

[\begin{matrix} - D & E_{1} & E_{2} \end{matrix}] [\begin{matrix} t \\ u \\ v \end{matrix}] = T

克莱姆法则¶

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

t = \frac{1}{| \begin{matrix} - D & E_{1} & E_{2} \end{matrix} |} | \begin{matrix} T & E_{1} & E_{2} \end{matrix} |

u = \frac{1}{| \begin{matrix} - D & E_{1} & E_{2} \end{matrix} |} | \begin{matrix} - D & T & E_{2} \end{matrix} |

v = \frac{1}{| \begin{matrix} - D & E_{1} & E_{2} \end{matrix} |} | \begin{matrix} - D & E_{1} & T \end{matrix} |

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

[\begin{matrix} t \\ u \\ v \end{matrix}] = \frac{1}{| \begin{matrix} - D & E_{1} & E_{2} \end{matrix} |} | \begin{matrix} T & E_{1} & E_{2} \\ - D & T & E_{2} \\ - D & E_{1} & T \end{matrix} |

向量混合积¶

根据混合积公式

| \begin{matrix} a & b & c \end{matrix} | = a \times b ∙ c = - a \times c ∙ b

上式可以改写成

[\begin{matrix} t \\ u \\ v \end{matrix}] = \frac{1}{[\begin{matrix} D \times E_{2} ∙ E_{1} \end{matrix}]} [\begin{matrix} T \times E_{1} ∙ E_{2} \\ D \times E_{2} ∙ T \\ T \times E_{1} ∙ D \end{matrix}]

令

P = D \times E_{2} Q = T \times E_{1}

最终的计算公式¶

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

[\begin{matrix} t \\ u \\ v \end{matrix}] = \frac{1}{| \begin{matrix} P ∙ E_{1} \end{matrix} |} | \begin{matrix} Q ∙ E_{2} \\ P ∙ T \\ Q ∙ D \end{matrix} |

代码¶

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;
}