雁过无痕

  C++博客 :: 首页 :: 新随笔 :: 联系 :: 聚合  :: 管理 ::

前一篇文章讨论了二维坐标系下的判断,下面讨论三维的情况。

三维坐标下,对点的判断明显要复杂很多。如果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 中国大陆许可协议进行许可,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接,否则保留追究法律责任的权利。
posted on 2011-07-14 23:28 flyinghearts 阅读(1469) 评论(0)  编辑 收藏 引用 所属分类: 算法

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理