计算机图形学——射线与三角形相交检测_Möller-Trumbore算法及推导过程(涉及标量三重积、克莱姆法则)
1. 标量三重积
1.1 标量三重积的定义
标量三重积(scalar triple product)是一个结合点积和叉积的运算,定义为:
a⋅(b×c) \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) a⋅(b×c)
这里:
- b×c\mathbf{b} \times \mathbf{c}b×c:是向量 b\mathbf{b}b 和 c\mathbf{c}c 的叉积,结果是一个垂直于 b\mathbf{b}b 和 c\mathbf{c}c 所在平面的向量。
- a⋅(b×c)\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c})a⋅(b×c):将 a\mathbf{a}a 投影到 b×c\mathbf{b} \times \mathbf{c}b×c 上后,计算标量乘积,最终结果是一个标量。
1.2 标量三重积的几何意义
1.2.1 平行六面体的体积
标量三重积的绝对值表示由三个向量 a,b,c\mathbf{a}, \mathbf{b}, \mathbf{c}a,b,c 作为边所构成的平行六面体的体积:
V=∣a⋅(b×c)∣ V = |\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c})| V=∣a⋅(b×c)∣
几何解释:

- b×c\mathbf{b} \times \mathbf{c}b×c 的模长代表 b\mathbf{b}b 和 c\mathbf{c}c 所围成的平行四边形的面积。
- a⋅(b×c)\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c})a⋅(b×c) 是平行四边形面积乘以向量 a\mathbf{a}a 在其法向上的投影,即体积。
1.2.2 向量共面的判别
若 a⋅(b×c)=0\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = 0a⋅(b×c)=0,说明向量 a,b,c\mathbf{a}, \mathbf{b}, \mathbf{c}a,b,c 共面。
1.2.3 方向判定
标量三重积的符号可用于判断 a,b,c\mathbf{a}, \mathbf{b}, \mathbf{c}a,b,c 的排列是否符合右手法则:
- 若 a⋅(b×c)>0\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) > 0a⋅(b×c)>0,则符合右手法则。
- 若 a⋅(b×c)<0\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) < 0a⋅(b×c)<0,则不符合。
1.3 标量三重积的计算公式
标量三重积可以通过行列式计算:
a⋅(b×c)=∣a1a2a3b1b2b3c1c2c3∣ \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = \begin{vmatrix} a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \\ c_1 & c_2 & c_3 \end{vmatrix} a⋅(b×c)= a1b1c1a2b2c2a3b3c3
推导过程:
- 叉积 b×c\mathbf{b} \times \mathbf{c}b×c 的结果为:
b×c=[ijkb1b2b3c1c2c3]=(b2c3−b3c2)i−(b1c3−b3c1)j+(b1c2−b2c1)k \mathbf{b} \times \mathbf{c} = \begin{bmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ b_1 & b_2 & b_3 \\ c_1 & c_2 & c_3 \end{bmatrix} = (b_2c_3 - b_3c_2)\mathbf{i} - (b_1c_3 - b_3c_1)\mathbf{j} + (b_1c_2 - b_2c_1)\mathbf{k} b×c= ib1c1jb2c2kb3c3 =(b2c3−b3c2)i−(b1c3−b3c1)j+(b1c2−b2c1)k
- 点积 a⋅(b×c)\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c})a⋅(b×c):
a⋅(b×c)=a1(b2c3−b3c2)+a2(b3c1−b1c3)+a3(b1c2−b2c1) \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = a_1(b_2c_3 - b_3c_2) + a_2(b_3c_1 - b_1c_3) + a_3(b_1c_2 - b_2c_1) a⋅(b×c)=a1(b2c3−b3c2)+a2(b3c1−b1c3)+a3(b1c2−b2c1)
这正是对应的 3×33 \times 33×3 行列式展开式。
1.4 标量三重积的性质
1.4.1 交换性
标量三重积具有循环交换的性质:
a⋅(b×c)=b⋅(c×a)=c⋅(a×b) \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = \mathbf{b} \cdot (\mathbf{c} \times \mathbf{a}) = \mathbf{c} \cdot (\mathbf{a} \times \mathbf{b}) a⋅(b×c)=b⋅(c×a)=c⋅(a×b)
但改变叉积的顺序会改变符号:
a⋅(b×c)=−a⋅(c×b)=−b⋅(a×c)=−c⋅(b×a) \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = -\mathbf{a} \cdot (\mathbf{c} \times \mathbf{b}) = -\mathbf{b} \cdot (\mathbf{a} \times \mathbf{c}) = -\mathbf{c} \cdot (\mathbf{b} \times \mathbf{a}) a⋅(b×c)=−a⋅(c×b)=−b⋅(a×c)=−c⋅(b×a)
1.4.2 零条件
若标量三重积为零,说明三个向量共面,因为平行六面体的体积为零。
1.4.3 分配律
标量三重积满足分配律:
a⋅[(b+c)×d]=a⋅(b×d)+a⋅(c×d) \mathbf{a} \cdot [(\mathbf{b} + \mathbf{c}) \times \mathbf{d}] = \mathbf{a} \cdot (\mathbf{b} \times \mathbf{d}) + \mathbf{a} \cdot (\mathbf{c} \times \mathbf{d}) a⋅[(b+c)×d]=a⋅(b×d)+a⋅(c×d)
1.5 标量三重积的应用
1.5.1 几何问题
- 平行六面体体积:利用标量三重积计算三维空间中任意三个向量张成的平行六面体体积。
- 向量共面性判定:通过标量三重积是否为零判断三个向量是否共面。
1.5.2 物理学
- 力矩:在力学中,标量三重积可以描述作用力和位置矢量与某方向的力矩关系。
- 电磁学:在电磁场分析中,标量三重积用于描述矢量场的空间关系。
1.5.3 线性代数
- 判断向量是否线性相关:若 a,b,c\mathbf{a}, \mathbf{b}, \mathbf{c}a,b,c 线性相关,则标量三重积为零。
- 计算行列式:三阶行列式的值与标量三重积一致。
2. 克莱姆法则(Cramer’s Rule)
克莱姆法则(Cramer’s Rule)是线性代数中的一种方法,用于求解具有唯一解的线性方程组,尤其是当系数矩阵为方阵时。它通过行列式来表示方程组的解。以下是详细的讲解。
2.1 克莱姆法则的背景和适用范围
克莱姆法则适用于以下线性方程组:
Ax=b A \mathbf{x} = \mathbf{b} Ax=b
其中:
- AAA 是 n×nn \times nn×n 的系数矩阵(方阵)。
- x\mathbf{x}x 是未知量的列向量 (x1,x2,…,xn)T(x_1, x_2, \dots, x_n)^T(x1,x2,…,xn)T。
- b\mathbf{b}b 是常数项的列向量 (b1,b2,…,bn)T(b_1, b_2, \dots, b_n)^T(b1,b2,…,bn)T。
前提条件:
- 矩阵 AAA 是方阵。
- 矩阵 AAA 的行列式 ∣A∣≠0|A| \neq 0∣A∣=0(保证方程组有唯一解)。
2.2 克莱姆法则的核心思想
克莱姆法则通过引入“替换矩阵”来求解每个未知数。公式如下:
xi=∣Ai∣∣A∣ x_i = \frac{|A_i|}{|A|} xi=∣A∣∣Ai∣
其中:
- ∣A∣|A|∣A∣ 是系数矩阵 AAA 的行列式。
- ∣Ai∣|A_i|∣Ai∣ 是将系数矩阵 AAA 的第 iii 列替换为常数向量 b\mathbf{b}b 后的新矩阵的行列式。
替换矩阵的构造:
矩阵 AiA_iAi 是通过将 AAA 的第 iii 列替换为 b\mathbf{b}b 得到的:
Ai=[a11…b1…a1na21…b2…a2n⋮…⋮…⋮an1…bn…ann] A_i = \begin{bmatrix} a_{11} & \dots & b_1 & \dots & a_{1n} \\ a_{21} & \dots & b_2 & \dots & a_{2n} \\ \vdots & \dots & \vdots & \dots & \vdots \\ a_{n1} & \dots & b_n & \dots & a_{nn} \end{bmatrix} Ai=
a11a21⋮an1…………b1b2⋮bn…………a1na2n⋮ann
2.3 克莱姆法则的推导
从矩阵方程 Ax=bA \mathbf{x} = \mathbf{b}Ax=b 开始:
- 解方程的核心是找到向量 x\mathbf{x}x,其表达式为:
x=A−1b \mathbf{x} = A^{-1} \mathbf{b} x=A−1b
其中 A−1A^{-1}A−1 是矩阵 AAA 的逆。
- 通过矩阵行列式的性质,逆矩阵 A−1A^{-1}A−1 可以用伴随矩阵表示:
A−1=adj(A)∣A∣ A^{-1} = \frac{\text{adj}(A)}{|A|} A−1=∣A∣adj(A)
其中 adj(A)\text{adj}(A)adj(A) 是 AAA 的伴随矩阵,∣A∣|A|∣A∣ 是 AAA 的行列式。
- 结合行列式计算,将伴随矩阵的行和 b\mathbf{b}b 相乘,得出每个 xi=∣Ai∣∣A∣x_i = \frac{|A_i|}{|A|}xi=∣A∣∣Ai∣ 的公式。
2.4 计算步骤
-
求系数矩阵 AAA 的行列式 ∣A∣|A|∣A∣:
如果 ∣A∣=0|A| = 0∣A∣=0,克莱姆法则不适用,因为方程组没有唯一解。 -
构造替换矩阵 AiA_iAi:
将 AAA 的第 iii 列用 b\mathbf{b}b 替换。 -
计算替换矩阵的行列式 ∣Ai∣|A_i|∣Ai∣。
-
求解未知量:
用公式 xi=∣Ai∣∣A∣x_i = \frac{|A_i|}{|A|}xi=∣A∣∣Ai∣ 计算每个 xix_ixi。
2.5 例子
考虑线性方程组:
{2x+y−z=8−3x−y+2z=−11−2x+y+2z=−3 \begin{cases} 2x + y - z = 8 \\ -3x - y + 2z = -11 \\ -2x + y + 2z = -3 \end{cases} ⎩ ⎨ ⎧2x+y−z=8−3x−y+2z=−11−2x+y+2z=−3
- 系数矩阵 AAA:
A=[21−1−3−12−212] A = \begin{bmatrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2 \end{bmatrix} A= 2−3−21−11−122
- 行列式 ∣A∣|A|∣A∣:
计算行列式:
∣A∣=∣21−1−3−12−212∣=2(−1⋅2−1⋅2)−1(−3⋅2−(−2)⋅2)+(−1)(−3⋅1−(−2)⋅(−1))=−6+16−5=5 |A| = \begin{vmatrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2 \end{vmatrix} = 2(-1 \cdot 2 - 1 \cdot 2) - 1(-3 \cdot 2 - (-2) \cdot 2) + (-1)(-3 \cdot 1 - (-2) \cdot (-1)) = -6 + 16 - 5 = 5 ∣A∣= 2−3−21−11−122 =2(−1⋅2−1⋅2)−1(−3⋅2−(−2)⋅2)+(−1)(−3⋅1−(−2)⋅(−1))=−6+16−5=5
- 构造替换矩阵 A1A_1A1, A2A_2A2, A3A_3A3:
- 替换第 1 列(求 x1x_1x1):
A1=[81−1−11−12−312] A_1 = \begin{bmatrix} 8 & 1 & -1 \\ -11 & -1 & 2 \\ -3 & 1 & 2 \end{bmatrix} A1= 8−11−31−11−122
- 替换第 2 列(求 x2x_2x2):
A2=[28−1−3−112−2−32] A_2 = \begin{bmatrix} 2 & 8 & -1 \\ -3 & -11 & 2 \\ -2 & -3 & 2 \end{bmatrix} A2= 2−3−28−11−3−122
- 替换第 3 列(求 x3x_3x3):
A3=[218−3−1−11−21−3] A_3 = \begin{bmatrix} 2 & 1 & 8 \\ -3 & -1 & -11 \\ -2 & 1 & -3 \end{bmatrix} A3= 2−3−21−118−11−3
-
计算替换矩阵的行列式:
- ∣A1∣=5|A_1| = 5∣A1∣=5,∣A2∣=10|A_2| = 10∣A2∣=10,∣A3∣=15|A_3| = 15∣A3∣=15。
-
求解解向量:
x1=∣A1∣∣A∣=55=1,x2=∣A2∣∣A∣=105=2,x3=∣A3∣∣A∣=155=3 x_1 = \frac{|A_1|}{|A|} = \frac{5}{5} = 1, \quad x_2 = \frac{|A_2|}{|A|} = \frac{10}{5} = 2, \quad x_3 = \frac{|A_3|}{|A|} = \frac{15}{5} = 3 x1=∣A∣∣A1∣=55=1,x2=∣A∣∣A2∣=510=2,x3=∣A∣∣A3∣=515=3
解得:
x=[123] \mathbf{x} = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix} x= 123
2.6 克莱姆法则的优点和局限性
优点:
- 公式化强,直接通过行列式计算解。
- 在手工计算小规模方程组时,方法清晰易懂。
局限性:
- 对于高维方程组(如 n>3n > 3n>3),行列式计算复杂,效率较低。
- 仅适用于系数矩阵为方阵且行列式不为零的情况,不能处理增广矩阵的特殊解或无解情况。
3. Möller–Trumbore 算法
Möller–Trumbore算法是一种高效的射线与三角形相交检测算法,由 Tomas Möller 和 Ben Trumbore 在1997年提出。这种算法在光线追踪(Ray Tracing)、碰撞检测(Collision Detection)以及计算机图形学中广泛使用,因其速度快、计算量小的特点而非常适合实时应用。
3.1 算法的主要目标
给定一条射线 R(t)=O+td\mathbf{R}(t) = \mathbf{O} + t\mathbf{d}R(t)=O+td 和一个三角形(顶点为 V0\mathbf{V}_0V0, V1\mathbf{V}_1V1, V2\mathbf{V}_2V2),判断射线是否与三角形相交。如果相交,返回:
- 交点的位置 ttt:表示交点沿射线的参数化距离。
- 交点的重心坐标 (1−b1−b2,b1,b2)(1 - b_1 - b_2, b_1, b_2)(1−b1−b2,b1,b2),用来判断交点是否在三角形内。
3.2 关键概念
-
三角形表示:
一个三角形由三个顶点 v0\mathbf{v}_0v0、v1\mathbf{v}_1v1、v2\mathbf{v}_2v2 定义。 -
光线表示:
光线由起点 O\mathbf{O}O 和方向向量 d\mathbf{d}d 定义:
R(t)=O+td,t≥0 \mathbf{R}(t) = \mathbf{O} + t\mathbf{d}, \quad t \geq 0 R(t)=O+td,t≥0
- 重心坐标:
Möller–Trumbore算法通过重心坐标插值的原理,将三角形上任意点P\mathbf{P}P 表示为:
P=(1−b1−b2)V0+b1V1+b2V2(1) \mathbf{P} = (1 - b_1 - b_2)V_0 + b_1V_1 + b_2V_2 \tag{1} P=(1−b1−b2)V0+b1V1+b2V2(1)
交点 P\mathbf{P}P 必须同时满足这两个方程,将射线与此表达式联立,转化为线性方程组并解出参数 ttt、b1b_1b1、b2b_2b2。
3.2 Möller-Trumbore算法推导过程
Step 1
设三角形的三个顶点为 V0,V1,V2V_0, V_1, V_2V0,V1,V2,根据重心坐标插值,我们知道三角形所在平面上任意一点与射线相交可以写成:
O+td=(1−b1−b2)V0+b1V1+b2V2(1) \mathbf{O}+ t\mathbf{d} = (1 - b_1 - b_2)V_0 + b_1V_1 + b_2V_2 \tag{1} O+td=(1−b1−b2)V0+b1V1+b2V2(1)
我们只要算出 b1,b2b_1, b_2b1,b2,根据重心坐标插值的规则,直接就可以知道交点在不在三角形内。
上面的式子只有三个未知数:t,b1,b2t, b_1, b_2t,b1,b2,下面的任务就是把它们解出来。
Step 2
将 (1) 式子变换为:
O−V0=(V1−V0)b1+(V2−V0)b2−td O - V_0 = (V_1 - V_0)b_1 + (V_2 - V_0)b_2 - t\mathbf{d} O−V0=(V1−V0)b1+(V2−V0)b2−td
定义:
S=O−V0,E1=V1−V0,E2=V2−V0 \mathbf{S} = O - V_0, \quad \mathbf{E}_1 = V_1 - V_0, \quad \mathbf{E}_2 = V_2 - V_0 S=O−V0,E1=V1−V0,E2=V2−V0
得到:
S=E1b1+E2b2−td \mathbf{S} = \mathbf{E}_1b_1 + \mathbf{E}_2b_2 - t\mathbf{d} S=E1b1+E2b2−td
将其写成矩阵的形式:
[−dE1E2][tb1b2]=S [-\mathbf{d} \quad \mathbf{E}_1 \quad \mathbf{E}_2] \begin{bmatrix} t \\ b_1 \\ b_2 \end{bmatrix} = \mathbf{S} [−dE1E2] tb1b2 =S
Step 3
根据克莱姆法则,解得 ttt:
t=det([SE1E2])det([−dE1E2]) t = \frac{\text{det}([\mathbf{S} \quad \mathbf{E}_1 \quad \mathbf{E}_2])}{\text{det}([-\mathbf{d} \quad \mathbf{E}_1 \quad \mathbf{E}_2])} t=det([−dE1E2])det([SE1E2])
根据三重积性质:
det([SE1E2])=(S×E1)⋅E2 \text{det}([\mathbf{S} \quad \mathbf{E}_1 \quad \mathbf{E}_2]) = (\mathbf{S} \times \mathbf{E}_1) \cdot \mathbf{E}_2 det([SE1E2])=(S×E1)⋅E2
det([−dE1E2])=−d⋅(E1×E2)=E1⋅(d×E2) \text{det}([-\mathbf{d} \quad \mathbf{E}_1 \quad \mathbf{E}_2]) = -\mathbf{d} \cdot (\mathbf{E}_1 \times \mathbf{E}_2) = \mathbf{E}_1 \cdot (\mathbf{d} \times \mathbf{E}_2) det([−dE1E2])=−d⋅(E1×E2)=E1⋅(d×E2)
因此得到:
t=(S×E1)⋅E2E1⋅(d×E2) t = \frac{(\mathbf{S} \times \mathbf{E}_1) \cdot \mathbf{E}_2}{\mathbf{E}_1 \cdot (\mathbf{d} \times \mathbf{E}_2)} t=E1⋅(d×E2)(S×E1)⋅E2
定义:
S1=d×E2,S2=S×E1 \mathbf{S}_1 = \mathbf{d} \times \mathbf{E}_2, \quad \mathbf{S}_2 = \mathbf{S} \times \mathbf{E}_1 S1=d×E2,S2=S×E1
可得:
t=S2⋅E2S1⋅E1 t = \frac{\mathbf{S}_2 \cdot \mathbf{E}_2}{\mathbf{S}_1 \cdot \mathbf{E}_1} t=S1⋅E1S2⋅E2
Step 4
同理可以得到 b1b_1b1 和 b2b_2b2:
b1=S1⋅SS1⋅E1,b2=S2⋅dS1⋅E1 b_1 = \frac{\mathbf{S}_1 \cdot \mathbf{S}}{\mathbf{S}_1 \cdot \mathbf{E}_1}, \quad b_2 = \frac{\mathbf{S}_2 \cdot \mathbf{d}}{\mathbf{S}_1 \cdot \mathbf{E}_1} b1=S1⋅E1S1⋅S,b2=S1⋅E1S2⋅d
整理总结:
其中:
E1=V1−V0\mathbf{E}_1 = \mathbf{V}_1 - \mathbf{V}_0E1=V1−V0
E2=V2−V0\mathbf{E}_2 = \mathbf{V}_2 - \mathbf{V}_0E2=V2−V0
S=O−V0\mathbf{S} = \mathbf{O} - \mathbf{V}_0S=O−V0
S1=D×E2\mathbf{S}_1 = \mathbf{D} \times\mathbf{E}_2S1=D×E2
S2=S×E1\mathbf{S}_2 = \mathbf{S} \times\mathbf{E}_1S2=S×E1
3.3 判断交点是否在三角形内
解得 t,b1,b2t, b_1, b_2t,b1,b2 后:
- 如果 t<0t < 0t<0,说明交点在射线起点的反方向。
- 如果 b1,b2≥0b_1, b_2 \geq 0b1,b2≥0 且 1−b1−b2≥01 - b_1 - b_2 \geq 01−b1−b2≥0,说明交点位于三角形内部。
3.4 光线与三角形的相交的代码实现
// 检测光线是否与三角形相交
// 参数:
// orig: 光线的起点(通常是相机位置或反射光线的起点)
// dir: 光线的方向(单位向量)
// v0, v1, v2: 三角形的三个顶点
// tnear: 输出参数,光线起点到交点的距离
// u, v: 输出参数,重心坐标中的两个参数,用于判断交点是否在三角形内部
// 返回值:如果光线与三角形相交,返回 true,并更新 tnear, u, v,否则返回 false
bool rayTriangleIntersect(const Vector3f &orig, const Vector3f &dir,
const Vector3f &v0, const Vector3f &v1, const Vector3f &v2,
float &tnear, float &u, float &v)
{
// 计算三角形的两条边
Vector3f E1 = v1 - v0;
Vector3f E2 = v2 - v0;
Vector3f S = orig - v0;
Vector3f S1 = crossProduct(dir, E2);
Vector3f S2 = crossProduct(S, E1);
// 计算行列式(determinant),用于判断光线是否平行于三角形
float det = dotProduct(S1, E1);
if (det <= 0) return false; // 如果行列式为 0 或小于 0,光线与三角形平行或背向三角形
u = dotProduct(S1, S) / det;
if (u < 0 || u > 1) return false; // 如果 u 不在 [0, 1] 范围内,交点在三角形外部
v = dotProduct(S2, dir) / det;
if (v < 0 || u + v > 1) return false; // 如果 v 不在 [0, 1 - u] 范围内,交点在三角形外部
// 计算光线参数 tnear
tnear = dotProduct(S2, E2) / det;
if (tnear < 0) return false; // 如果 tnear 小于 0,交点在光线的反方向
return true; // 光线与三角形相交,返回 true
}
3.5 Möller-Trumbore算法的应用场景
- 光线追踪:判断光线是否击中物体表面。
- 碰撞检测:快速检测点是否落在多边形表面。
- 三角网格渲染:确定屏幕上像素是否位于三角形内。
DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。
更多推荐


所有评论(0)