1. 概述
任何复杂的三维模型都可以视作空间三角面片的集合,很容易碰到的一个问题就是空间射线与三角形相交的问题,例如拾取、遮蔽检测等。这里就总结下该问题的两种算法实现。
2. 常规算法
一种很常规的思路就是先计算射线与三角面片的交点,再看该交点是否再三角形内部。
2.1. 理论推导
对于空间一条射线,令起点为O,其方向为D,根据射线的参数公式,其上任意一点P(也就是要求的交点)为:
P
=
O
+
t
D
(1)
P = O + tD \tag{1}
P=O+tD(1)
其中t>0,根据t的取值不同,可得射线上不同的点,也就是关键在于求未知量t的值。
已知空间三角面片三个顶点为v1,v2,v3,那么很容易可以求得三角面片的法向量n。显然面上的向量(v1-P)与n是垂直的,则它们的点积为0:
(
v
1
−
P
)
⋅
n
=
0
(2)
(v1-P) \cdot n = 0 \tag{2}
(v1−P)⋅n=0(2)
将式(1)代入式(2),求得未知量t为:
t
=
(
v
1
−
O
)
⋅
n
D
⋅
n
t = \frac{ (v1-O) \cdot n }{ D \cdot n}
t=D⋅n(v1−O)⋅n
再将t代入到(1)式中,即可得到射线与该三点组成的平面了。
接下来就是判断这个交点是否在三角形面之内了,由于是空间三角形,所以比较好的算法是文献[2]中提到的同向法,摘录如下:
![同向法](https://i-blog.csdnimg.cn/blog_migrate/9a01ab7a8d7139154af03b22684ee3ea.png)
2.2. 具体实现
具体的C/C++实现代码如下:
#include <iostream>
using namespace std;
#define EPSILON 0.000001
// 3D vector
class Vector3d
{
public:
Vector3d()
{
}
~Vector3d()
{
}
Vector3d(double dx, double dy, double dz)
{
x = dx;
y = dy;
z = dz;
}
// 矢量赋值
void set(double dx, double dy, double dz)
{
x = dx;
y = dy;
z = dz;
}
// 矢量相加
Vector3d operator + (const Vector3d& v) const
{
return Vector3d(x + v.x, y + v.y, z + v.z);
}
// 矢量相减
Vector3d operator - (const Vector3d& v) const
{
return Vector3d(x - v.x, y - v.y, z - v.z);
}
//矢量数乘
Vector3d Scalar(double c) const
{
return Vector3d(c*x, c*y, c*z);
}
// 矢量点积
double Dot(const Vector3d& v) const
{
return x * v.x + y * v.y + z * v.z;
}
// 矢量叉积
Vector3d Cross(const Vector3d& v) const
{
return Vector3d(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
}
double _x()
{
return x;
}
double _y()
{
return y;
}
double _z()
{
return z;
}
private:
double x, y, z;
};
// v1 = Cross(AB, AC)
// v2 = Cross(AB, AP)
// 判断矢量v1和v2是否同向
bool SameSide(Vector3d A, Vector3d B, Vector3d C, Vector3d P)
{
Vector3d AB = B - A;
Vector3d AC = C - A;
Vector3d AP = P - A;
Vector3d v1 = AB.Cross(AC);
Vector3d v2 = AB.Cross(AP);
// v1 and v2 should point to the same direction
//return v1.Dot(v2) >= 0 ;
return v1.Dot(v2) > 0;
}
// 判断点P是否在三角形ABC内(同向法)
bool PointinTriangle1(Vector3d A, Vector3d B, Vector3d C, Vector3d P)
{
return SameSide(A, B, C, P) && SameSide(B, C, A, P) && SameSide(C, A, B, P);
}
//ray-triangle intersection algorithm (通过平面方程计算)
//参数说明:V1,V2,V3,三角形三点;O,射线原点;D,射线方向
bool ray_triangle_intersection1(Vector3d V1, Vector3d V2, Vector3d V3, Vector3d O, Vector3d D, Vector3d *I)
{
bool rv = false;
//v1(n1,n2,n3);
//平面方程: na * (x – n1) + nb * (y – n2) + nc * (z – n3) = 0 ;
double na = (V2._y() - V1._y())*(V3._z() - V1._z()) - (V2._z() - V1._z())*(V3._y() - V1._y());
double nb = (V2._z() - V1._z())*(V3._x() - V1._x()) - (V2._x() - V1._x())*(V3._z() - V1._z());
double nc = (V2._x() - V1._x())*(V3._y() - V1._y()) - (V2._y() - V1._y())*(V3._x() - V1._x());
//平面法向量
Vector3d nv(na, nb, nc);
//平面法向量与射线方向向量差积
double vpt = D.Dot(nv);
if (vpt == 0)
{
rv = false; //此时直线与平面平行
}
else
{
Vector3d P = V1 - O;
double t = P.Dot(nv) / vpt;
*I = O + D.Scalar(t);
if (PointinTriangle1(V1, V2, V3, *I))
{
rv = true;
}
else
{
rv = false;
}
}
return rv;
}
int main()
{
Vector3d V1(0, 0, 0);
Vector3d V2(50, 0, 0);
Vector3d V3(0, 50, 0);
Vector3d O(5, 10, -10);
Vector3d P(10, 10, 10);
Vector3d D = P - O;
Vector3d I;
if (ray_triangle_intersection1(V1, V2, V3, O, D, &I)) {
cout << I._x() << '\t' << I._y() << '\t' << I._z() << endl;
}
}
3. 优化算法
仔细思考常规算法的思路,在计算射线与平面的交点的时候,实际是将射线的参数方程与平面的参数方程联立求值即可。那么如果知道空间三角形的参数方程,将其与射线的参数方程联立,不就可以直接求得交点了吗?Tomas Moller的论文《Fast, Minimum Storage Ray Triangle Intersection》提出了一种优化算法,正是基于这个思路,并且给出了合理的解法。
3.1. 理论推导
对于三个顶点为V1,V2,V3组成的空间三角形,对于三角形内的任一点,有如下参数方程:
P
=
(
1
−
u
−
v
)
V
1
+
u
V
2
+
v
V
3
(3)
P = (1-u-v)V1 + uV2 + vV3 \tag{3}
P=(1−u−v)V1+uV2+vV3(3)
u, v是V2和V3的权重,1-u-v是V1的权重,并且满足u>=0, v >= 0,u+v<=1。这个参数方程的具体解释可参考文献[5],摘录如下:
![三角形参数方程](https://i-blog.csdnimg.cn/blog_migrate/1c18a19dceec298c10076a7a810d8027.png)
将射线公式(1)与三角形公式(3)联立起来,有:
(
1
−
u
−
v
)
V
1
+
u
V
2
+
v
V
3
=
O
+
t
D
(1-u-v)V1 + uV2 + vV3 = O + tD
(1−u−v)V1+uV2+vV3=O+tD
很显然,u、v、t都是未知数,移项并整理,可得如下线性方程组:
[
−
D
V
2
−
V
1
V
3
−
V
1
]
[
t
u
v
]
=
O
−
V
1
\left[ \begin{matrix} -D & V2-V1 & V3-V1 \end{matrix} \right] \left[ \begin{matrix} t \\ u \\ v \end{matrix} \right] = O - V1
[−DV2−V1V3−V1]⎣⎡tuv⎦⎤=O−V1
可以使用克莱姆法则来求解这个线性方程组,大家可以复习下线性代数(文献[6]),我这里也将其摘录如下:
![克莱姆法则](https://i-blog.csdnimg.cn/blog_migrate/9d1284d95d783be5e1510ee87d095a6b.png)
令
E
1
=
V
2
−
V
1
,
E
2
=
V
3
−
V
1
,
T
=
O
−
V
1
E1 = V2 - V1,E2 = V3 - V1,T = O - V1
E1=V2−V1,E2=V3−V1,T=O−V1,则上式可以改写成:
[
−
D
E
1
E
2
]
[
t
u
v
]
=
T
\left[ \begin{matrix} -D & E1 & E2 \end{matrix} \right] \left[ \begin{matrix} t \\ u \\ v \end{matrix} \right] = T
[−DE1E2]⎣⎡tuv⎦⎤=T
根据克莱姆法则,有:
{
t
=
1
∣
−
D
E
1
E
2
∣
∣
T
E
1
E
2
∣
u
=
1
∣
−
D
E
1
E
2
∣
∣
−
D
T
E
2
∣
v
=
1
∣
−
D
E
1
E
2
∣
∣
−
D
E
1
T
∣
\begin{cases} t = \frac{1}{\begin{vmatrix} -D&E1&E2\\ \end{vmatrix}} \begin{vmatrix} T&E1&E2\\ \end{vmatrix} \\ u = \frac{1}{\begin{vmatrix} -D&E1&E2\\ \end{vmatrix}} \begin{vmatrix} -D&T&E2\\ \end{vmatrix} \\ v = \frac{1}{\begin{vmatrix} -D&E1&E2\\ \end{vmatrix}} \begin{vmatrix} -D&E1&T\\ \end{vmatrix} \\ \end{cases}
⎩⎪⎪⎪⎨⎪⎪⎪⎧t=∣−DE1E2∣1∣∣TE1E2∣∣u=∣−DE1E2∣1∣∣−DTE2∣∣v=∣−DE1E2∣1∣∣−DE1T∣∣
接下来就要用到向量的混合积公式(具体可参看文献[7])了,对于三向量a,b,c,有:
∣
a
b
c
∣
=
a
×
b
⋅
c
=
−
a
×
c
⋅
b
=
−
c
×
b
⋅
a
\begin{vmatrix} a&b&c\\ \end{vmatrix} = a \times b \cdot c = - a \times c \cdot b = - c \times b \cdot a
∣∣abc∣∣=a×b⋅c=−a×c⋅b=−c×b⋅a
上式可改写成:
{
t
=
1
D
×
E
2
⋅
E
1
(
T
×
E
1
⋅
E
2
)
u
=
1
D
×
E
2
⋅
E
1
(
D
×
E
2
⋅
T
)
v
=
1
D
×
E
2
⋅
E
1
(
T
×
E
1
⋅
D
)
\begin{cases} t = \frac{1}{D \times E2 \cdot E1} (T \times E1 \cdot E2) \\ u = \frac{1}{D \times E2 \cdot E1} (D \times E2 \cdot T) \\ v = \frac{1}{D \times E2 \cdot E1} (T \times E1 \cdot D) \\ \end{cases}
⎩⎪⎨⎪⎧t=D×E2⋅E11(T×E1⋅E2)u=D×E2⋅E11(D×E2⋅T)v=D×E2⋅E11(T×E1⋅D)
令
P
=
D
×
E
2
,
Q
=
T
×
E
1
P=D \times E2, Q = T \times E1
P=D×E2,Q=T×E1,进一步简化可得:
{
t
=
1
P
⋅
E
1
(
Q
⋅
E
2
)
u
=
1
P
⋅
E
1
(
P
⋅
T
)
v
=
1
P
⋅
E
1
(
Q
⋅
D
)
\begin{cases} t = \frac{1}{P \cdot E1} (Q \cdot E2) \\ u = \frac{1}{P \cdot E1} (P \cdot T) \\ v = \frac{1}{P \cdot E1} (Q \cdot D) \\ \end{cases}
⎩⎪⎨⎪⎧t=P⋅E11(Q⋅E2)u=P⋅E11(P⋅T)v=P⋅E11(Q⋅D)
3.2. 具体实现
具体的C/C++实现代码如下:
#include <iostream>
using namespace std;
#define EPSILON 0.000001
// 3D vector
class Vector3d
{
public:
Vector3d()
{
}
~Vector3d()
{
}
Vector3d(double dx, double dy, double dz)
{
x = dx;
y = dy;
z = dz;
}
// 矢量赋值
void set(double dx, double dy, double dz)
{
x = dx;
y = dy;
z = dz;
}
// 矢量相加
Vector3d operator + (const Vector3d& v) const
{
return Vector3d(x + v.x, y + v.y, z + v.z);
}
// 矢量相减
Vector3d operator - (const Vector3d& v) const
{
return Vector3d(x - v.x, y - v.y, z - v.z);
}
//矢量数乘
Vector3d Scalar(double c) const
{
return Vector3d(c*x, c*y, c*z);
}
// 矢量点积
double Dot(const Vector3d& v) const
{
return x * v.x + y * v.y + z * v.z;
}
// 矢量叉积
Vector3d Cross(const Vector3d& v) const
{
return Vector3d(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
}
double _x()
{
return x;
}
double _y()
{
return y;
}
double _z()
{
return z;
}
private:
double x, y, z;
};
//ray-triangle intersection algorithm
//参数说明:V1,V2,V3,三角形三点;O,射线原点;D,射线方向。
bool ray_triangle_intersection(Vector3d V1, Vector3d V2, Vector3d V3, Vector3d O, Vector3d D, Vector3d *I)
{
//Find vectors for two edges sharing V1
Vector3d e1 = V2 - V1;
Vector3d e2 = V3 - V1;
//Begin calculating determinant - also used to calculate u parameter
Vector3d P = D.Cross(e2);
//if determinant is near zero, ray lies in plane of triangle
double det = e1.Dot(P);
//NOT CULLING
if (det > -EPSILON && det < EPSILON)
{
return false;
}
double inv_det = 1.f / det;
//calculate distance from V1 to ray origin
Vector3d T = O - V1;
//Calculate u parameter and test bound
double u = T.Dot(P) * inv_det;
//The intersection lies outside of the triangle
if (u < 0.f || u > 1.f)
{
return false;
}
//Prepare to test v parameter
Vector3d Q = T.Cross(e1);
//Calculate V parameter and test bound
double v = D.Dot(Q) * inv_det;
//The intersection lies outside of the triangle
if (v < 0.f || u + v > 1.f)
{
return false;
}
double t = e2.Dot(Q) * inv_det;
//ray intersection
if (t > EPSILON)
{
*I = O + D.Scalar(t);
return true;
}
return false;
}
int main()
{
Vector3d V1(0, 0, 0);
Vector3d V2(50, 0, 0);
Vector3d V3(0, 50, 0);
Vector3d O(5, 10, -10);
Vector3d P(10, 10, 10);
Vector3d D = P - O;
Vector3d I;
if (ray_triangle_intersection(V1, V2, V3, O, D, &I)) {
cout << I._x() << '\t' << I._y() << '\t' << I._z() << endl;
}
}
可以看到这种优化算法无论是代码量还是时间、空间复杂度都由于原来的常规算法,最直观的体现就是判断语句多,能够即使返回避免后续运算。
4. 参考
[1] Möller–Trumbore intersection algorithm
[2] 判断点是否在三角形内
[3] 射线与平面的相交检测(Ray-Plane intersection test)
[4] 射线和三角形的相交检测(ray triangle intersection test)
[5] 三角形方程? - 高崎汀步的回答 - 知乎
[6] 克莱姆法则
[7] 三矢量的混合积