空間射線與三角形相交演算法的兩種實現
目錄
- 1. 概述
- 2. 常規演算法
- 2.1. 理論推導
- 2.2. 具體實現
- 3. 優化演算法
- 3.1. 理論推導
- 3.2. 具體實現
- 4. 參考
1. 概述
任何複雜的三維模型都可以視作空間三角面片的集合，很容易碰到的一個問題就是空間射線與三角形相交的問題，例如拾取、遮蔽檢測等。這裡就總結下該問題的兩種演算法實現。
2. 常規演算法
一種很常規的思路就是先計算射線與三角面片的交點，再看該交點是否再三角形內部。
2.1. 理論推導
對於空間一條射線，令起點為O，其方向為D，根據射線的引數公式，其上任意一點P（也就是要求的交點）為：
\[P = O + tD \tag{1}\]
其中t>0，根據t的取值不同，可得射線上不同的點，也就是關鍵在於求未知量t的值。
已知空間三角面片三個頂點為v1，v2，v3，那麼很容易可以求得三角面片的法向量n。顯然面上的向量(v1-P)與n是垂直的，則它們的點積為0：
\[(v1-P) \cdot n = 0 \tag{2}\]
將式(1)代入式(2)，求得未知量t為：
\[t = \frac{ (v1-O) \cdot n }{ D \cdot n} \]
再將t代入到(1)式中，即可得到射線與該三點組成的平面了。
接下來就是判斷這個交點是否在三角形面之內了，由於是空間三角形，所以比較好的演算法是文獻[2]中提到的同向法，摘錄如下：
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)V1 + uV2 + vV3 \tag{3}\]
u, v是V2和V3的權重，1-u-v是V1的權重，並且滿足u>=0, v >= 0,u+v<=1。這個引數方程的具體解釋可參考文獻[5]，摘錄如下：
將射線公式(1)與三角形公式(3)聯立起來，有：
\[(1-u-v)V1 + uV2 + vV3 = O + tD\]
很顯然，u、v、t都是未知數，移項並整理，可得如下線性方程組：
\[
\left[
\begin{matrix}
-D & V2-V1 & V3-V1
\end{matrix}
\right]
\left[
\begin{matrix}
t \\
u \\
v
\end{matrix}
\right] = O - V1
\]
可以使用克萊姆法則來求解這個線性方程組，大家可以複習下線性代數(文獻[6])，我這裡也將其摘錄如下：
令\(E1 = V2 - V1，E2 = V3 - V1，T = O - V1\)，則上式可以改寫成：
\[
\left[
\begin{matrix}
-D & E1 & E2
\end{matrix}
\right]
\left[
\begin{matrix}
t \\
u \\
v
\end{matrix}
\right] = 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}
\]
接下來就要用到向量的混合積公式（具體可參看文獻[7]）了，對於三向量a，b，c，有：
\[
\begin{vmatrix}
a&b&c\\
\end{vmatrix} = a \times b \cdot c = - a \times c \cdot b = - c \times b \cdot a
\]
上式可改寫成：
\[\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}
\]
令\(P=D \times E2, Q = T \times E1\)，進一步簡化可得：
\[\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}
\]
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] 三向量的混