设为首页 加入收藏

TOP

点在三角形内之二 (三维坐标系1)
2014-11-23 22:30:33 来源: 作者: 【 】 浏览:1
Tags:三角形 之二 三维 坐标系

三维坐标下,对点的判断明显要复杂很多。如果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 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 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次))。

】【打印繁体】【投稿】【收藏】 【推荐】【举报】【评论】 【关闭】 【返回顶部
分享到: 
上一篇HDU 1301 Jungle Roads 下一篇HDU 1233 还是畅通工程

评论

帐  号: 密码: (新用户注册)
验 证 码:
表  情:
内  容: