点在三角形ABC内 之二 (三维坐标系1).
前一篇文章讨论了二维坐标系下的判断,下面讨论三维的情况。
三维坐标下,对点的判断明显要复杂很多。如果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次))。
作者: flyinghearts
出处: http://www.cnblogs.com/flyinghearts/
本文采用知识共享署名-非商业性使用-相同方式共享 2.5 中国大陆许可协议进行许可,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接,否则保留追究法律责任的权利。