前一篇文章討論了二維坐標系下的判斷，下面討論三維的情況。

三維坐標下，對點的判斷明顯要復(fù)雜很多。如果google“Point in triangle”，第一個搜索結(jié)果就是這個網(wǎng)頁，可惜的是，作者沒有對結(jié)果進一步討論，沒有給出一個好的實現(xiàn)，而且其所有結(jié)論只在已知四點共面時才成立。

 

 

 

 

前面已經(jīng)證明過，面積法和向量同向法是等價的。

ab × ac = ab × ap + ap × ac + pb × pc

|ab × ac| = |ab × ap| + |ap × ac| + |pb × pc|

由于ab × ac、ab × ap、ap × ac、pb × pc這4個向量平行同向，因而可以先判斷a、b、c、p這四點是否共面（通過計算混合積），若共面的話，則4個向量一定共線，接著只要判斷后三個向量是否都和第一個向量（ab × ac）同向（可以通過判斷后三個向量與第一個向量的點積的正負性來確定）。

 

代碼：

 

template<typename T> class Vec3 {

 T x, y, z;

public:

 Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {}

 Vec3 operator-(const Vec3& v) const { return Vec3(x - v.x, y - v.y, z - v.z); }

 T dot(const Vec3& v) const { return x * v.x + y * v.y + z * v.z; }

 

 Vec3 cross(const Vec3& v) const {

    return Vec3(y * v.z - z * v.y, z * v.x - x * v.z,x * v.y - y * v.x);

 }

};

 

typedef Vec3<double> V3d;

 

//方法一

bool is_in_triangle3a(const V3d& a, const V3d& b, const V3d& c, const V3d& p)

{

 V3d ab(b - a), ac(c - a), ap(p -a);

 if (fabs(ab.cross(ac).dot(ap)) >= DBL_EPSILON) return false; //四點不共面

 V3d abc = ab.cross(ac), abp = ab.cross(ap), apc = ap.cross(ac);

 double t0 = abc.dot(abc), t1 = abp.dot(abc), t2 = apc.dot(abc);

 

 //t1 >= 0     t2 >= 0    t1 + t2 <= t0       t0肯定大于0

 // return (t1 >= -DBL_EPSILON) & (t2 >= -DBL_EPSILON) & (t0 - t1 - t2 >= -DBL_EPSILON);

 double delta = fabs(t1) + fabs(t2) + fabs(t0 - t1 - t2) - t0;

  return fabs(delta) < DBL_EPSILON; 

}

 

方法一，需要30次乘法計算。即使在已知四點共面的情況下，仍需要27次乘法計算，僅節(jié)省了3次乘法計算?？紤]到每次計算向量積需要6次乘法計算，而計算點積只要3次乘法計算，因而可以考慮消除向量積計算：

 

 

利用公式：

 (a × b)· c = (c × a)· b　　　    （混合積）

 a × (b × c) = b（a·c）− c（a·b）  (拉格朗日公式)

可得：

(a × b)·(a × c) = ((a × c) × a)·b = ((a·a)c – (a·c)a)·b

= (a·a) * (b·c) – (a·c) * (a·b)

 

利用這個展開式，可得：

 

//方法二

bool is_in_triangle3b(const V3d& a, const V3d& b, const V3d& c, const V3d& p)

{

 V3d ab(b - a), ac(c - a), ap(p -a);

 if (fabs(ab.cross(ac).dot(ap)) >= DBL_EPSILON) return false; //四點不共面

 V3d abc = ab.cross(ac), abp = ab.cross(ap), apc = ap.cross(ac);

 

 //double t0 = abc.dot(abc), t1 = abp.dot(abc), t2 = apc.dot(abc); //對這三個計算公式進行展開

 double v11 = ab.dot(ab), v22 = ac.dot(ac), v12 = ab.dot(ac);

 double v13 = ab.dot(ap), v23 = ac.dot(ap);

 double t0 = v11 * v22 - v12 * v12;

 double t1 = v11 * v23 - v12 * v13;

 double t2 = v22 * v13 - v12 * v23;

 

 double delta = fabs(t1) + fabs(t2) + fabs(t0 - t1 - t2) - t0;

 return fabs(delta) < DBL_EPSILON; 

}

 

方法二，需要30次乘法計算，但在已知四點共面時則只需要21次乘法計算。

 

上面的兩種方法，方法一，容易記，容易實現(xiàn)，且在不能確定四點共面時，效率與方法二差不多（甚至可能略高）；而方法二最大的優(yōu)點，則是在已知四點共面時，比方法一少用6次乘法，但是實現(xiàn)起來實在麻煩。那么，是否存在更好的方法呢？答案是肯定的，這就是后面要提到的方法（多一次條件判斷，只要13次乘法（四點共面時只要8次））。

 

作者： flyinghearts
出處： http://www.cnblogs.com/flyinghearts/
本文采用知識共享署名-非商業(yè)性使用-相同方式共享 2.5 中國大陸許可協(xié)議進行許可，歡迎轉(zhuǎn)載，但未經(jīng)作者同意必須保留此段聲明，且在文章頁面明顯位置給出原文連接，否則保留追究法律責(zé)任的權(quán)利。