前一篇文章讨论了二维坐标系下的判断，下面讨论三维的情况。
三维坐标下，对点的判断明显要复杂很多。如果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
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次））。
最新评论