射线与三角形的交点

前言

射线与三角形的求交是一道经典面试题，它是模型选择的必备基础。

这道题要是说的话，并不难，但要往深了说，也不太容易。

而你若想脱颖而出，那就得往深了说。

我接下来会整体的说一下这个解题的思路和方法，有些公式我不会说太细，大家哪里不懂可以点击后面的链接学习。

前期准备

两种解法

射线与三角形求交的常见解法有两种：

先求交，后过滤：先求射线与三角形所在平面的交点，然后判断此交点是否在三角形中。
先过滤，后求交：逐步过滤射线与三角形相交的条件，最后推导出交点。

第一种方法比较简单，很容易推导和理解，但计算效率会比较低，需要具备的知识是叉乘、点积、射线与平面的求交公式。

第二种方法相对难一些，需要具备扎实的线性变换基础，计算效率会高于前者，需要具备的知识是叉乘、点积、矩阵行列式、克莱姆法则(Cramer’s rule)。

咱们依次说一下这两种方法。

解法一

我们用”先求交，后过滤”的方法求射线与三角形交点。

已知：

三角形ABC
射线，原点为O，方向为d

求：射线与三角形的交点P

解：

1.通过AB和AC的叉乘，求出三角形ABC所在平面的垂线n。

若n等于0，则三角形的三点共线，甚至共点，需要过滤掉。

注：这里的n 不需要归一化。

2.通过射线与平面的求交公式求出交点。

根据射线与平面的求交公式可得：

t：射线上的时间或距离标量
O：射线原点
A：平面上一点

射线的方程是：

r(t) ：在射线上到射线原点的距离为t的点位。

将之前求出的t值代入r(t)便可求出交点，设此交点为P。

3.判断点P是否在三角形ABC中。

在右手坐标系中，若AP^AC,BP^BA,CP^CB分别点积AB^AC都大于0，则点P在三角形ABC中，此时的点P就是射线与三角形的交点，如下图所示：

否则，点P在三角形ABC之外，如下图所示：

解法二

我们用”先过滤，后求交”的方法求射线与三角形交点。

以具象化的方式把问题和问题的解法画出来，要比数字计算更容易理解。所以我接下来会把这个问题的解法画出来。

这个解法的核心是：基于射线的方向和三角形的两条边把射线的源点投影到一个新的空间中。

已知：

三角形ABC
射线，原点为O，方向为d

求：射线与三角形的交点P

解：

根据已知条件，我们可以得到一个线性变换矩阵[-d,AB,AC]，矩阵也可以理解为坐标系。

在这个矩阵中，A是原点，-d是射线的反方向，AB,AC是三角形的两条边。

向量AO在坐标系[-d,AB,AC] 中的本地坐标位是(t,u,v)。

由上面的条件可知，(t,u,v) 就是向量AO在[-d,AB,AC] 中的本地坐标位。

所以(t,u,v) 就是[-d,AB,AC]的逆矩阵乘以AO。

因为-d是一个单位向量，所以t值就是世界坐标系中射线与三角形ABC相交的距离，由此距离可以得到交点P：

现在我们已经用一个很快的方式把射线与三角形的交点P给算出来了，那说好的过滤呢？

过滤就是要在刚才的计算过程中进行逐步拦截，判断射线有没有可能和三角形存在交点，若不存在，就不用再往后走了。

一边过滤一边求交

1.先判断矩阵[-d,AB,AC]的有效性。

[-d,AB,AC]的有效性可以通过其行列式来判断。

三维矩阵的行列式可以理解为由以矩阵基向量为临边的平行六面体的有向体积。

[-d,AB,AC] 的行列式的计算方式如下：

det 是determinant 的缩写，即行列式
AB叉乘AC的结果是垂直于AB、AC的向量，此向量的长度是以AB、AC为临边的平行四边形的面积
-d点积上面的向量，就是把-d投影到AB和AC的垂线上，从而得到平行六面体的高，然后乘以上面的向量的向量的长度，也就是平行四边形的面积，从而就得到了以-d,AB,AC为临边的平行六面体的有向体积，效果如下：

det([-d,AB,AC]) 的值有三种情况：

为正时，射线从三角形所在平面的正面穿过；
为负时，射线从三角形所在平面的背面穿过，若三角形背面不可选，则无交点；
为零时，矩阵[-d,AB,AC]发生空间降维，射线可能与三角形平行或者出现零向量，射线与三角形所在的平面没有交点或者有无数个交点。此情况一般不考虑，默认无交点。

2.把[-d,AB,AC]的逆矩阵乘以AO的过程分解、分别计算u,v,t，并过滤。

这时候就需要用到克莱姆法则，其规则如下：

上面的结构和咱们之前求(t,u,v) 时的结构是一样的：

顺便回顾一下我们在解法一里说过的，射线与平面的求交公式：

大家有没有觉得这两种解法里的t 有点相似？对，它们不仅相似，而且一样。

对于为什么用克莱姆法则就能算出t,u,v 来，我会在最后简单解释一下。

当然大家也可以看我在前期准备里提供的链接。

接下来咱们先往后走，说一下过滤条件。

当t 满足以下条件时不会有交点：

t 是从射线的原点O 向着射线的方向d 推进的距离，当此值小于0时，是向射线反方向推进的，这肯定是不对的。

当u和v 满足以下任一条件都不会有交点：

u<0
v<0
u+v>1

(t,u,v)可以理解为在[-d,AB,AC] 中的本地坐标。

当u或v小于0 时，交点会跑到三角形外面，这个比较好理解，因为点A就是源点，[-d,AB,AC] 中的基向量AB、AC 是三角形的边。

那u+v≤1 是什么概念呢？

这就是一个简单的向量加法。

想象你在直角边为1的等腰直角三角形内跑步，三角形中的任一点都可以理解为你先沿着AB轴跑了一段，然后又沿着AC轴跑了一段，你能跑出的最远距离只能是直角边的长度

3.把t 代入射线方程求出交点P。

代码实现

这是我用three.js写的一个计算过程。

php复制代码/* 三角形ABC */
const A=new Vector3(0, 0, 0,)
const B=new Vector3(0, 0, -1)
const C=new Vector3(1, 0, -1)
// 背面是否可选
let backfaceCulling=false

/* 射线 */
// 射线原点
const O=new Vector3(0,1,0)
//射线方向
const d=new Vector3(0.2,-1,-0.8)

/* 射线与三角形的求交 */
// 射线的反方向
const _d=new Vector3(-d.x,-d.y,-d.z)
// 向量AB
const AB=new Vector3().subVectors( B, A )
// 向量AC
const AC=new Vector3().subVectors( C, A )
// AB和AC垂线，其长度是以AB,AC为临边的平行四边形的面积
const n=new Vector3().crossVectors( AB, AC )
// 矩阵[-d,AB,AC]的行列式，即以-d,AB,AC为临边的平行六面体的有向体积
let det = _d.dot( n )
if ( det < 0 ) {
  if ( backfaceCulling ){
    return null
  }
} else if ( det === 0 ) {
  return null;
}
// 点A到射线源点的向量
const AO=new Vector3().subVectors( O, A )
// 从射线的原点O 向着射线的方向d 推进的距离
const t=AO.dot( n )/det
// 当t<0的时候，向射线反方向推进，与实际不符
if ( t < 0 ) {
  return null
}
// AO在AB上的投影坐标
const u =  _d.dot(new Vector3().crossVectors( AO, AC ))/det;
if ( u < 0 ) {
  return null
}
// AO在AC上的投影坐标
const v =  _d.dot(new Vector3().crossVectors( AB, AO ) )/det;
if ( v < 0 ) {
  return null
}
// 当u + v >1时会超出三角形的范围
if ( u + v >1  ) {
  return null
}
// r ( t ) = O + t d
const p=O.clone().add(d.clone().multiplyScalar(t))

在three.js的Ray对象里有一个intersectTriangle() 方法。

这个方法略有冗余，比如其中的sign 代表了行列式的取向，但这是不需要的，因为我们可以在两个行列式相除的时候，得到行列式取向负负得正的结果。

kotlin复制代码intersectTriangle( a, b, c, backfaceCulling, target ) {
        // Compute the offset origin, edges, and normal.
        // from https://github.com/pmjoniak/GeometricTools/blob/master/GTEngine/Include/Mathematics/GteIntrRay3Triangle3.h
        _edge1.subVectors( b, a );
        _edge2.subVectors( c, a );
        _normal.crossVectors( _edge1, _edge2 );

        // Solve Q + t*D = b1*E1 + b2*E2 (Q = kDiff, D = ray direction,
        // E1 = kEdge1, E2 = kEdge2, N = Cross(E1,E2)) by
        //   |Dot(D,N)|*b1 = sign(Dot(D,N))*Dot(D,Cross(Q,E2))
        //   |Dot(D,N)|*b2 = sign(Dot(D,N))*Dot(D,Cross(E1,Q))
        //   |Dot(D,N)|*t = -sign(Dot(D,N))*Dot(Q,N)
        let DdN = this.direction.dot( _normal );
        let sign;
        if ( DdN > 0 ) {
            if ( backfaceCulling ) return null;
            sign = 1;
        } else if ( DdN < 0 ) {
            sign = - 1;
            DdN = - DdN;
        } else {
            return null;
        }
        _diff.subVectors( this.origin, a );
        const DdQxE2 = sign * this.direction.dot( _edge2.crossVectors( _diff, _edge2 ) );

        // b1 < 0, no intersection
        if ( DdQxE2 < 0 ) {
            return null;
        }

        const DdE1xQ = sign * this.direction.dot( _edge1.cross( _diff ) );
        // b2 < 0, no intersection
        if ( DdE1xQ < 0 ) {
            return null;
        }

        // b1+b2 > 1, no intersection
        if ( DdQxE2 + DdE1xQ > DdN ) {
            return null;
        }

        // Line intersects triangle, check if ray does.
        const QdN = - sign * _diff.dot( _normal );

        // t < 0, no intersection
        if ( QdN < 0 ) {
            return null;
        }

        // Ray intersects triangle.
        return this.at( QdN / DdN, target );
    }