空间中直线段和三角形的相交算法

最近在看recast&detour源码的时候有遇到许多数学上的算法问题，特此记录，以便以后查看。

分析

问题： 已经线段pq (p起点 q终点) ，三角形abc（a b c逆时针存储），判断pq与abc有无交点。

第一步：

判断三角形所在的面和线段所在的直线 是否平行 或者 是否是从三角形面的背面进入。如果是就提前退出，否则进入第二步。

（使用面的法向量 和 线段的向量 的 点积）

第二步：

判断线段和三角形所在的面的交点是否在线段上。

平面的方程可以表示为：

     X*normal = d                                              （1）

其中 X 为平面上的点，d为原点到平面的距离。

以p为起点的，沿着pq方向的射线上的点的坐标可以表示为：

M = p + t*pq = p - t*qp                                    （2）

假设M点即为 平面与pq所在的直线的交点：(M点和a点均为在平面上的点，带入公式（1）中

\[(\overrightarrow {op}  - t \cdot \overrightarrow {qp} ) \cdot \overrightarrow {normal}  = \overrightarrow {oa}  \cdot \overrightarrow {normal} \]

解得t为

\[t = \frac{{\overrightarrow {ap}  \cdot \overrightarrow {normal} }}{{\overrightarrow {qp}  \cdot \overrightarrow {normal} }}\]

判断t是否 （0,1）之间，如果是继续第三步，如果否返回false。

第三步：

判断是否交点在三角形内部。

三角形abc内的点M可以用三角形的重心坐标系表示 ：

\[M = (1 - {\lambda _2} - {\lambda _3})a + {\lambda _2}b + {\lambda _3}c\]

其中如果lamda3 和 lamda2 都在（0,1）之间，那么表示M在abc的内部

改写成向量形式，且把M用公式（2）带入，得到：

\[\overrightarrow {op}  - t\overrightarrow {qp}  = (1 - {\lambda _2} - {\lambda _3})\overrightarrow {oa}  + {\lambda _2}\overrightarrow {ob}  + {\lambda _3}\overrightarrow {oc} \]

整理得：

\[{\lambda _2}\overrightarrow {ab}  + {\lambda _3}\overrightarrow {ac}  + t\overrightarrow {qp}  = \overrightarrow {ap} \]

写成向量的形式：

\[\left( {\begin{array}{*{20}{c}}
  {a{b_x}}&{a{c_x}}&{q{p_x}} \\ 
  {a{b_y}}&{a{c_y}}&{q{p_y}} \\ 
  {a{b_z}}&{a{c_z}}&{q{p_z}} 
\end{array}} \right)\left( {\begin{array}{*{20}{c}}
  {{\lambda _2}} \\ 
  {{\lambda _3}} \\ 
  t 
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
  {a{p_x}} \\ 
  {a{p_y}} \\ 
  {a{p_z}} 
\end{array}} \right)\]

求解lamda2 lamda3：（应用tips中的克莱姆法则以及混合积的运算法则）

\[\begin{gathered}
  {\lambda _2} = \frac{{(ap,ac,qp)}}{{(ab,ac,qp)}} = \frac{{ac \cdot (qp \times ap)}}{{qp \cdot (ab \times ac)}} \hfill \\
  {\lambda _3} = \frac{{(ab,ap,qp)}}{{(ab,ac,qp)}} = \frac{{ab \cdot (ap \times qp)}}{{qp \cdot (ab \times ac)}} = \frac{{ - ab \cdot (qp \times ap)}}{{qp \cdot (ab \times ac)}} \hfill \\ 
\end{gathered} \]

 

源码

// 空间点 sp 起点 sq终点
// 三角形空间点 a b c 
// 输出参数 t 
static bool intersectSegmentTriangle(const float* sp, const float* sq,
									 const float* a, const float* b, const float* c,
									 float &t)
{
	float v, w;
	float ab[3], ac[3], qp[3], ap[3], norm[3], e[3];
	rcVsub(ab, b, a);
	rcVsub(ac, c, a);
	// 终点->起点 的向量, 不是 起点->终点
	rcVsub(qp, sp, sq);
	
	// Compute triangle normal. Can be precalculated or cached if
	// intersecting multiple segments against the same triangle
	// 法线方向为三角形正面方向 
	rcVcross(norm, ab, ac);
	
	// Compute denominator d. If d <= 0, segment is parallel to or points
	// away from triangle, so exit early
	// d = 0.0 说明 qp 和 norm 垂直，说明三角形和 qp 平行。
	// d < 0.0 说明 qp 和 norm 是钝角 说明是从三角形的背面 进入和三角形相交的 
	float d = rcVdot(qp, norm);
	if (d <= 0.0f) return false;

	// Compute intersection t value of pq with plane of triangle. A ray
	// intersects if 0 <= t. Segment intersects if 0 <= t <= 1. Delay
	// dividing by d until intersection has been found to pierce triangle
	rcVsub(ap, sp, a);
	t = rcVdot(ap, norm);
	if (t < 0.0f) return false;
	if (t > d) return false; // For segment; exclude this code line for a ray test 
	
	// Compute barycentric coordinate components and test if within bounds
	rcVcross(e, qp, ap);
	v = rcVdot(ac, e);
	if (v < 0.0f || v > d) return false;
	w = -rcVdot(ab, e);
	if (w < 0.0f || v + w > d) return false;
	
	// Segment/ray intersects triangle. Perform delayed division
	t /= d;
	
	return true;
}

tips

克莱姆法则

一个线性方程组可以用矩阵（A为n*n的方阵）与向量（x，c均长度为n的列向量）的方程来表示：
\[Ax = c\]

如果  A是一个可逆矩阵，那么方程有解，为

\[{x_i} = \frac{{\det ({A_i})}}{{\det (A)}}\]

其中Ai是被列向量 c取代了 A的第i列的列向量后得到的矩阵。

混合积的运算法则

1.混合积和三维矩阵行列式的关系（其中 后一个等号 矩阵的转置矩阵的行列式等于这个矩阵的行列式。）

\[a \cdot (b \times c) = \left| {\begin{array}{*{20}{c}}
  {{a_1}}&{{a_2}}&{{a_3}} \\ 
  {{b_1}}&{{b_2}}&{{b_3}} \\ 
  {{c_1}}&{{c_2}}&{{c_3}} 
\end{array}} \right| = \left| {\begin{array}{*{20}{c}}
  {{a_1}}&{{b_1}}&{{c_1}} \\ 
  {{a_2}}&{{b_2}}&{{c_2}} \\ 
  {{a_3}}&{{b_3}}&{{c_3}} 
\end{array}} \right|\]

2.混合积的交换律

\[a \cdot (b \times c) = b \cdot (c \times a) = c \cdot (a \times b)\]