点三角的测试
1.同方技术
判断三维空间中的点是否在一个三角形中,最简单的方法就是进行向量叉积运算,怎么,没懂?那就接着往下看,假设空间中由这样一个三角形,黄色区域代表三角形内部,如何判断空间中的任意一点是否在黄色区域中呢?
或许你还没有想到比较简单快速的方法,再来看看下面的这张图,或许你会豁然开朗:
以上图为例,取三角形外的一点p,向量p-A与向量B-A的叉积的方向指向屏幕外;取三角形内的一点p’,向量p’-A与向量B-A的叉积方向指向屏幕,由此我们得出直线AB上方的点的矢量,与矢量B-A的叉积指向屏幕外,而直线AB下方的点的矢量,得到的叉积指向屏幕,这样我们利用叉积就能确定平面上的点位于直线AB的哪一侧。但是一般情况下,我们如何确定叉积应该指向那个方向呢,这时,我们发现点C还没用到,如果点p’再三角形的内部,那么p’与C应该位于直线AB的同侧,也就是两个点分别与向量B-A的叉积大方向是同一个方向。同样的方法,如果我们判断出点p’与A位于直线BC的同一侧,点p’与B位于直线AC的同一侧,那么就能保证点p’位于黄色区域内部,即位于三角形内部。代码实现也很简单:
function SameSide(p1,p2,a,b)
cp1 = CrossProduct(ba,p1-a)
cp2 = CrossProduct(ba,p2-a)
if DotProduct(cp1,cp2)> = 0
return true
else return false
function PointInTriangle( p,a,b,c)
if SameSide(p,a,b,c)and SameSide(p,b,a,c ) and SameSide(p,c,a,b)
return true
else return false
对向量叉积知识不熟悉的同学,我在这里补充一下:
两向量a,b的叉积表示为axb,但要跟数学中的乘法区分开来!其表达式为:
两向量叉积axb同时垂直a,b:
这种重心方法简单有效,没有其他多余的运算。
2.重心技术
上面的方法核心就是确保点需要位于线的同一侧即可,它确实相当简单了。但是这就足够了吗,下面的重心技术不仅简单,而且执行速度更快!
我们都知道三角形的三个点在空间中定义了一个平面,选择其中任意一个点,我们可以将平面上所有其他位置视为相对于该点。若选择点A为起点,以向量(C-A)和向量(B-A)作为基础向量,就可以表示平面上的任意一点,只需要从点A开始沿着(C-A)走一段距离,再沿着(B-A)走一段距离即可。
这样,平面上的任意一点的位置就可以描述为:
P = A + u *(C - A)+ v *(B - A)
这里需要注意的是,u和v应该大于0,当u <0或v <0时,表示我们走错了方向,这个点位于三角形之外;如果u >1或v >1,表示我们沿着一个方向走向了三角形之外;最后,如果u+v>1,我们再次越过了边缘BC离开三角形。
给定u和v,我们就可以很轻松的用上面的方程计算出点P,但是如何能根据给定的点P计算出u和v呢,我们接着往下看:
P = A + u *(C - A)+ v *(B - A)//原方程
(P - A)= u *(C - A)+ v *(B - A)//两者边同时减去A
v2 = u * v0 + v * v1 //替换v0,v1,v2以减少写入
//我们有两个未知数(u和v)所以我们需要两个方程来解它们。在方程两边同时点乘v0得到第一个方程,然后在方程两边同时点乘v1得到第二个方程
(v2) . v0 = (u * v0 + v * v1) . v0
(v2) . v1 = (u * v0 + v * v1) . v1
//分配v0和v1
v2 . v0 = u * (v0 . v0) + v * (v1 . v0)
v2 . v1 = u * (v0 . v1) + v * (v1 . v1)
//现在我们有两个方程和两个未知数,可以解出一个变量并替换为另一个变量
求解 [v2.v0 == {u(v0.v0) + v(v1.v0), v2.v1 == u(v0.v1) + v(v1.v1)}, {u, v}]
u = ((v1.v1)(v2.v0)-(v1.v0)(v2.v1)) / ((v0.v0)(v1.v1) - (v0.v1)(v1.v0))
v = ((v0.v0)(v2.v1)-(v0.v1)(v2.v0)) / ((v0.v0)(v1.v1) - (v0.v1)(v1.v0))
算法表示为:
// Compute vectors
v0 = C - A
v1 = B - A
v2 = P - A
// Compute dot products
dot00 = dot(v0, v0)
dot01 = dot(v0, v1)
dot02 = dot(v0, v2)
dot11 = dot(v1, v1)
dot12 = dot(v1, v2)
// Compute barycentric coordinates
invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
u = (dot11 * dot02 - dot01 * dot12) * invDenom
v = (dot00 * dot12 - dot01 * dot02) * invDenom
// Check if point is in triangle
return (u >= 0) && (v >= 0) && (u + v < 1)