文章目錄

• 1. 概述
• 2. 常規算法
• • 2.1. 理論推導
  • 2.2. 具體實現
• 3. 優化算法
• • 3.1. 理論推導
  • 3.2. 具體實現
• 4. 參考

1. 概述

任何復雜的三維模型都可以視作空間三角面片的集合，很容易碰到的一個問題就是空間射線與三角形相交的問題，例如拾取、遮蔽檢測等。這里就總結下該問題的兩種算法實現。

2. 常規算法

一種很常規的思路就是先計算射線與三角面片的交點，再看該交點是否再三角形內部。

2.1. 理論推導

對于空間一條射線，令起點為O，其方向為D，根據射線的參數公式，其上任意一點P（也就是要求的交點）為：
$$P=O+tD\text{\hspace{1em}(1)}$$

其中t>0，根據t的取值不同，可得射線上不同的點，也就是關鍵在于求未知量t的值。

已知空間三角面片三個頂點為v1，v2，v3，那么很容易可以求得三角面片的法向量n。顯然面上的向量(v1-P)與n是垂直的，則它們的點積為0：
$$(v1?P)?n=0\text{\hspace{1em}(2)}$$

將式(1)代入式(2)，求得未知量t為：
$$t=\frac{(v1?O)?n}{D?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\text{\hspace{1em}(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{array}{ccc}?D& V2?V1& V3?V1\end{array}\right]?\begin{array}{c}t\\ u\\ v\end{array}?=O?V1$$

可以使用克萊姆法則來求解這個線性方程組，大家可以復習下線性代數(文獻[6])，我這里也將其摘錄如下：

令$E1=V2?V1，E2=V3?V1，T=O?V1$，則上式可以改寫成：
$$\left[\begin{array}{ccc}?D& E1& E2\end{array}\right]?\begin{array}{c}t\\ u\\ v\end{array}?=T$$

根據克萊姆法則，有：
$$?\begin{array}{l}t=\frac{1}{|\begin{array}{ccc}?D& E1& E2\end{array}|}\left|\begin{array}{ccc}T& E1& E2\end{array}\right|\\ u=\frac{1}{|\begin{array}{ccc}?D& E1& E2\end{array}|}\left|\begin{array}{ccc}?D& T& E2\end{array}\right|\\ v=\frac{1}{|\begin{array}{ccc}?D& E1& E2\end{array}|}\left|\begin{array}{ccc}?D& E1& T\end{array}\right|\end{array}$$

接下來就要用到向量的混合積公式（具體可參看文獻[7]）了，對于三向量a，b，c，有：
$$\left|\begin{array}{ccc}a& b& c\end{array}\right|=a\times b?c=?a\times c?b=?c\times b?a$$
上式可改寫成：
$$?\begin{array}{l}t=\frac{1}{D\times E2?E1}(T\times E1?E2)\\ u=\frac{1}{D\times E2?E1}(D\times E2?T)\\ v=\frac{1}{D\times E2?E1}(T\times E1?D)\end{array}$$
令$P=D\times E2,Q=T\times E1$，進一步簡化可得：
$$?\begin{array}{l}t=\frac{1}{P?E1}(Q?E2)\\ u=\frac{1}{P?E1}(P?T)\\ v=\frac{1}{P?E1}(Q?D)\end{array}$$

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 V1Vector3d e1 = V2 - V1;Vector3d e2 = V3 - V1;//Begin calculating determinant - also used to calculate u parameterVector3d P = D.Cross(e2);//if determinant is near zero, ray lies in plane of triangledouble det = e1.Dot(P);//NOT CULLINGif (det > -EPSILON && det < EPSILON){return false;}double inv_det = 1.f / det;//calculate distance from V1 to ray originVector3d T = O - V1;//Calculate u parameter and test bounddouble u = T.Dot(P) * inv_det;//The intersection lies outside of the triangleif (u < 0.f || u > 1.f){return false;}//Prepare to test v parameterVector3d Q = T.Cross(e1);//Calculate V parameter and test bounddouble v = D.Dot(Q) * inv_det;//The intersection lies outside of the triangleif (v < 0.f || u + v > 1.f){return false;}double t = e2.Dot(Q) * inv_det;//ray intersectionif (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] 三矢量的混合積

總結

以上是生活随笔為你收集整理的空间射线与三角形相交算法的两种实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。