直線-三角形相交演算法總結

演算法一：直接與平面相交 + 交點是否處於三角形內部

直線與三角形相交，首先三角形是在平面上，直線與平面的交點比較容易求得。

平面公式：

直線可以寫成：

代入：

當我們計算出來 t 之後，再代回原公式得到交點P：

這裡有一點點危險的地方可能為0，這個時候就意味著 N 和 AB 垂直，那麼直線與平面平行，它與平面沒有交點。

當直線與平面相交之後，我們需要判斷交點是否在三角形內，之前在二維平面上，我們用過這樣的辦法：判斷和三角形邊的法向量的 n 的點乘是否小於0來看P點是否位於三角形邊的『左側』。

實際上這個思路也可以推廣到三維：

n 是三角形的normal，也是平面的normal是結果。如果P點在三角形內側，那麼的向量 n 會方向相同，也就是它們的點乘會大於0.

查看另外兩邊，也會有類似的效果，所以三角形和射線相交可以這樣寫：

bool ray_triangle_intersect(const Vec3f A, const Vec3f B, const Vec3f p0, const Vec3f p1, const Vec3f p2){

// 1. find the p point
Vec3f n = cross(p1 - p0, p2 - p0);
float d = -n*p1;

Vec3f AB = B - A;
if (fabsf(n*AB) < epsilon) return false;

float t = (-d - n*A )/(n * AB);
// t < 0 means opposite direction of AB
if (t < 0 ) return false;
Vec3f p = A + AB*t;

if (cross(p1-p0, p-p0) * n < 0 ) return false;
if (cross(p2-p1, p-p1) * n < 0 ) return false;
if (cross(p0-p2,p-p2) * n < 0 ) return false;
return true;
}

當然我們交點P也是直接求出來了。

演算法二：重心坐標系

重心坐標系我之前也寫過 → 三角形重心坐標

實際上重心坐標系也叫面積坐標，如果三角形面積為1的話，那麼P點分割的三角形有以下性質：

如果想驗證也很簡單：

所以這給我們提供了一個比較簡單的重心坐標系的演算法，再結合演算法一,我們可以很容易的算出 u, v.

// AB: ray, p0p1p2: triangle bool bary_centric_coordinate(const Vec3f A, const Vec3f B, const Vec3f p0, const Vec3f p1, const Vec3f p2, float &u, float &v){

// 1. find the p point, same as above
Vec3f n = cross(p1 - p0, p2 - p0);
float area = n.norm();
float d = -n*p1;

Vec3f AB = B - A;
if (fabsf(n*AB) < epsilon) return false;

float t = (-d - n*A )/(n * AB);
if (t < 0 ) return false;
Vec3f p = A + AB*t;

// 2. find u,v
Vec3f uvector, vvector;
if (cross(p2-p1, p-p1) * n < 0 ) return false;
if ((vvector = cross(p1-p0,p-p0)) * n < 0) return false;
if ((uvector = cross(p0-p2,p-p2)) * n < 0) return false;

u = uvector.norm()/area;
v = vvector.norm()/area;

return true;
}

實際上代碼跟演算法一很多部分都是一致的，只是這裡我們加了計算 u, v 而已。

M?ller-Trumbore 演算法

M?ller-Trumbore演算法我感覺和重心坐標系差不多，只是加入了更多線性代數的優化？三角形依舊是 ABC，假設光線源點為O，方向為D，有相交點滿足：

矩陣形式：

利用Cramers rule，可知：

其中 T = O - A, E1 = B - A, E2 = C - A.

又線性代數中行列式的性質：

繼續：

其中 .

所以整個計算就是我們無需再去計算三角形平面的一些性質，取而代之我們用以上式子就可以計算出 t, u, v.

// OD: ray, p0p1p2: triangle bool ray_triangle_intersect_mt(const Vec3f O, const Vec3f D, const Vec3f p0, const Vec3f p1, const Vec3f p2, float &t, float &u, float &v){ Vec3f e1 = p1-p0; Vec3f e2 = p2-p0; Vec3f pvec = cross(D, e2); float det = e1*pvec; if (det < epsilon) return false;