给定三角形ABC和一点P(x,y,z),判断点P是否在ABC内。这是游戏设计中一个常见的问题。需要注意的是,这里假定点和三角形位于同一个平面内。
内角和法
连接点P和三角形的三个顶点得到三条线段PA,PB和PC,求出这三条线段与三角形各边的夹角,如果所有夹角之和为360度,那么点P在三角形内,否则不在,此法直观,但效率低下。
同向法
假设点P位于三角形内,会有这样一个规律,当我们沿着ABCA的方向在三条边上行走时,你会发现点P始终位于边AB,BC和CA的右侧。我们就利用这一点,但是如何判断一个点在线段的左侧还是右侧呢?我们可以从另一个角度来思考,当选定线段AB时,点C位于AB的右侧,同理选定BC时,点A位于BC的右侧,最后选定CA时,点B位于CA的右侧,所以当选择某一条边时,我们只需验证点P与该边所对的点在同一侧即可。问题又来了,如何判断两个点在某条线段的同一侧呢?可以通过叉积来实现,连接PA,将PA和AB做叉积,再将CA和AB做叉积,如果两个叉积的结果方向一致,那么两个点在同一测。判断两个向量的是否同向可以用点积实现,如果点积大于0,则两向量夹角是锐角,否则是钝角。
代码如下,为了实现程序功能,添加了一个Vector3类,该类表示三维空间中的一个向量。
01 |
// 3D vector |
02 |
class Vector3
|
03 |
{ |
04 |
public :
|
05 |
Vector3( float fx, float fy, float fz)
|
06 |
:x(fx), y(fy), z(fz)
|
07 |
{
|
08 |
09 |
}
|
10 |
11 |
|
12 |
// Subtract
|
13 |
Vector3 operator - ( const Vector3& v) const
|
14 |
{
|
15 |
return Vector3(x - v.x, y - v.y, z - v.z) ;
|
16 |
}
|
17 |
18 |
|
19 |
// Dot product
|
20 |
float Dot( const Vector3& v) const
|
21 |
{
|
22 |
return x * v.x + y * v.y + z * v.z ;
|
23 |
}
|
24 |
25 |
|
26 |
// Cross product
|
27 |
Vector3 Cross( const Vector3& v) const
|
28 |
{
|
29 |
return Vector3(
|
30 |
y * v.z - z * v.y,
|
31 |
z * v.x - x * v.z,
|
32 |
x * v.y - y * v.x ) ;
|
33 |
}
|
34 |
35 |
36 |
public :
|
37 |
float x, y, z ;
|
38 |
}; |
39 |
40 |
41 |
// Determine whether two vectors v1 and v2 point to the same direction |
42 |
// v1 = Cross(AB, AC) |
43 |
// v2 = Cross(AB, AP) |
44 |
bool SameSide(Vector3 A, Vector3 B, Vector3 C, Vector3 P)
|
45 |
{ |
46 |
Vector3 AB = B - A ;
|
47 |
Vector3 AC = C - A ;
|
48 |
Vector3 AP = P - A ;
|
49 |
50 |
51 |
Vector3 v1 = AB.Cross(AC) ;
|
52 |
Vector3 v2 = AB.Cross(AP) ;
|
53 |
54 |
55 |
56 |
// v1 and v2 should point to the same direction
|
57 |
return v1.Dot(v2) >= 0 ;
|
58 |
} |
59 |
60 |
61 |
// Same side method |
62 |
// Determine whether point P in triangle ABC |
63 |
bool PointinTriangle1(Vector3 A, Vector3 B, Vector3 C, Vector3 P)
|
64 |
{ |
65 |
return SameSide(A, B, C, P) &&
|
66 |
SameSide(B, C, A, P) &&
|
67 |
SameSide(C, A, B, P) ;
|
68 |
} |
重心法
上面这个方法简单易懂,速度也快,下面这个方法速度更快,只是稍微多了一点数学而已
三角形的三个点在同一个平面上,如果选中其中一个点,其他两个点不过是相对该点的位移而已,比如选择点A作为起点,那么点B相当于在AB方向移动一段距离得到,而点C相当于在AC方向移动一段距离得到。
所以对于平面内任意一点,都可以由如下方程来表示:
P = A + u * (C – A) + v * (B - A) // 方程1
如果系数u或v为负值,那么相当于朝相反的方向移动,即BA或CA方向。那么如果想让P位于三角形ABC内部,u和v必须满足什么条件呢?有如下三个条件:
1 |
u >= 0 |
2 |
v >= 0 |
3 |
u + v <= 1 |
几个边界情况,当u = 0且v = 0时,就是点A,当u = 0,v = 1时,就是点B,而当u = 1, v = 0时,就是点C。
整理方程1得到P – A = u(C - A) + v(B - A)。
令v0 = C – A, v1 = B – A, v2 = P – A,则v2 = u * v0 + v * v1,现在是一个方程,两个未知数,无法解出u和v,将等式两边分别点乘v0和v1的到两个等式。
1 |
(v2) • v0 = (u * v0 + v * v1) • v0 |
2 |
(v2) • v1 = (u * v0 + v * v1) • v1 |
注意到这里u和v是数,而v0,v1和v2是向量,所以可以将点积展开得到下面的式子。
1 |
v2 • v0 = u * (v0 • v0) + v * (v1 • v0) // 式1
|
2 |
v2 • v1 = u * (v0 • v1) + v * (v1• v1) // 式2
|
解这个方程得到:
1 |
u = ((v1•v1)(v2•v0)-(v1•v0)(v2•v1)) / ((v0•v0)(v1•v1) - (v0•v1)(v1•v0)) |
2 |
v = ((v0•v0)(v2•v1)-(v0•v1)(v2•v0)) / ((v0•v0)(v1•v1) - (v0•v1)(v1•v0)) |
是时候上代码了,这段代码同样用到上面的Vector3类:
01 |
// Determine whether point P in triangle ABC |
02 |
bool PointinTriangle(Vector3 A, Vector3 B, Vector3 C, Vector3 P)
|
03 |
{ |
04 |
Vector3 v0 = C - A ;
|
05 |
Vector3 v1 = B - A ;
|
06 |
Vector3 v2 = P - A ;
|
07 |
08 |
|
09 |
float dot00 = v0.Dot(v0) ;
|
10 |
float dot01 = v0.Dot(v1) ;
|
11 |
float dot02 = v0.Dot(v2) ;
|
12 |
float dot11 = v1.Dot(v1) ;
|
13 |
float dot12 = v1.Dot(v2) ;
|
14 |
15 |
|
16 |
float inverDeno = 1 / (dot00 * dot11 - dot01 * dot01) ;
|
17 |
18 |
|
19 |
float u = (dot11 * dot02 - dot01 * dot12) * inverDeno ;
|
20 |
if (u < 0 || u > 1) // if u out of range, return directly
|
21 |
{
|
22 |
return false ;
|
23 |
}
|
24 |
25 |
|
26 |
float v = (dot00 * dot12 - dot01 * dot02) * inverDeno ;
|
27 |
if (v < 0 || v > 1) // if v out of range, return directly
|
28 |
{
|
29 |
return false ;
|
30 |
}
|
31 |
32 |
|
33 |
return u + v <= 1 ;
|
34 |
关于空间三角形的重心的解释: 在三角形ΔABC中,P点可以用ABC三点来表示,即P=sA+qB+rC 且s+q+r=1。实际上在由ABC三点所定义的平面上的任意点P都可以用ABC三点来表示,当0<s<1 0<q<1 0<r<1时点P落在三角形内部,即射线与三角形相交。 求解点P。 先前看过一个解法,来自MIT Open Course的Ray Tracing(Fall 2003)课程,大致方法如下: 因为 s+q+r=1,故有s=1-q-r,由此得P(q,r) = (1-q-r)A+qB+rC,化简得P(q,r) = A–(B-A)q+(C-A)r,现假设眼睛(射线起点)位置坐标围R,观察(射线)方向为D,则此射线上任意一点可以表示为P(t)=R+tD,假如三角形上的一点与射线相交,则有P(t)=P(q,r),即R+tD=A-(B-A)q+(C-A)r,将三个坐标分量带进去可得线性方程组: Rx+tDx=Ax-(Bx-Ax)q+(Cx-Ax)r Ry+tDy=Ay-(By-Ay)q+(Cy-Ay)r Rz+tDz=Az-(Bz-Az)q+(Cz-Az)r 写成矩阵形式即 由克莱姆法则即可求解三个变量 当q+r<1 且q>0r>0时交点即在三角形内部。 |A|=-[(A-B)ⅹ(C-A)]·D 当|A|>10-5 And |A|<105时返回 但是在实践中我发现一个问题,就是由s=1-q-r求出的s存在精度不高的问题,由克莱姆求解过程中涉及太多的乘法导致数据精度损失过大的问题,最终实际上s+q+r≠1,大约是s+q+r=1.000005,所以我们应该由三个方程变量分别为s,q,r构成的方程组来求出精度会比较高,并且尽量减少乘法的次数以保证精度。 用面积来求s,q,r: 由上面的知识我们知道s=Ts/T,q=Tq/T,r=Tr/T,而由给定射线及三角形三顶点可以直接求出焦点坐标。 1.射线上任意一点P(t) = R+tD; 2.假如P(t)在三角形平面上,则应满足平面方程(PX-X0)A+(PY-Y0)B+(PZ-Z0)C=0; 3.其中该平面的法线N=(NA,NB,NC),(X0,Y0,Z0)是平面上任意点; 4.将P(t)=R+tD带进方程即可得 (NARX+NBRY+NCRZ)+(NAtDX+NBtDY+NCtDZ)-(NACX+NBCY+NCCZ)=0; 5.或者表示为N·R+N·tD-N·C=0; 6.所t=(N·C-N·R)/N·D,法线N可由(A-B)ⅹ(A-C)求得,此方法并不会执行背面剔出因为法线方向不确定。 7.因为A×B=|A||B|sinθ=SΔABC/2 8.所以Ts/2=Normalize((B-P)ⅹ(C-P)); 9.s=Ts/Normalize(AⅹB) 9.用类似的方法还可以单独求出Tq,Tr 因为是后来才看到《实时计算机图形学》上关于射线三角形相交的算法介绍,相比较于后面自己弄出来的面积算法效率要高不少,乘法运算次数更少,在求解t的时候我已经消耗掉3个点乘,但是还是作为自己思考的一个结果写出来吧。这个方法不仅在光线追踪中很重要而且在实时游戏中也扮演重要角色,用途很广。 |