三維坐標下,對點的判斷明顯要復雜很多。如果google“Point in triangle”,第一個搜索結果就是這個網頁,可惜的是,作者沒有對結果進一步討論,沒有給出一個好的實現,而且其所有結論只在已知四點共面時才成立。
前面已經證明過,面積法和向量同向法是等價的。
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)同向(可以通過判斷後三個向量與第一個向量的點積的正負性來確定)。
代碼:
#include<cfloat>
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次乘法計算,僅節省了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次乘法計算。
上面的兩種方法,方法一,容易記,容易實現,且在不能確定四點共面時,效率與方法二差不多(甚至可能略高);而方法二最大的優點,則是在已知四點共面時,比方法一少用6次乘法,但是實現起來實在麻煩。那麼,是否存在更好的方法呢?答案是肯定的,這就是後面要提到的方法(多一次條件判斷,只要13次乘法(四點共面時只要8次))。