用矢量的叉积判断直线段是否有交
用矢量的叉积判断直线段是否有交
矢量叉积计算的另一个常用用途是直线段求交。求交算法是计算机图形学的核心算法,也是体现速度和稳定性的重要标志,高效并且稳定的求交算法是任何一个 CAD软件都必需要重点关注的。求交包含两层概念,一个是判断是否相交,另一个是求出交点。直线(段)的求交算法相对来说是比较简单的,首先来看看如何判 断两直线段是否相交。
常规的代数计算通常分三步,首先根据线段还原出两条线段所在直线的方程,然后联立方程组求出交点,最后再判断交点是否在线段区间上。常规的代数方法非常 繁琐,每次都要解方程组求交点,特别是交点不在线段区间的情况,计算交点就是做无用功。计算几何方法判断直线段是否有交点通常分两个步骤完成,这两个步骤 分别是快速排斥试验和跨立试验。假设要判断线段P1P2和线段Q1Q2是否有交点,则:
(1) 快速排斥试验
设以线段 P1P2 为对角线的矩形为R1, 设以线段 Q1Q2 为对角线的矩形为R2,如果R1和R2不相交,则两线段不会有交点;
(2) 跨立试验。
如 果两线段相交,则两线段必然相互跨立对方,所谓跨立,指的是一条线段的两端点分别位于另一线段所在直线的两边。判断是否跨立,还是要用到矢量叉积的几何意 义。以图3为例,若P1P2跨立Q1Q2 ,则矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的两侧,即:
( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0
上式可改写成:
( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0
当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时,说明线段P1P2和Q1Q2共线(但是不一定有交点)。同理判断Q1Q2跨立P1P2的依据是:
( Q1 - P1 ) × ( P2 - P1 ) * ( Q2 - P1 ) × ( P2 - P1 ) < 0
具体情况如下图所示:
图3 直线段跨立试验示意图
根据矢量叉积的几何意义,跨立试验只能证明线段的两端点位于另一个线段所在直线的两边,但是不能保证是在另一直线段的两端,因此,跨立试验只是证明两条 线段有交点的必要条件,必需和快速排斥试验一起才能组成直线段相交的充分必要条件。根据以上分析,两条线段有交点的完整判断依据就是:1)以两条线段为对 角线的两个矩形有交集;2)两条线段相互跨立。
判断直线段跨立用计算叉积算法的CrossProduct()函数即可,还需要一个判断两个矩形是否有交的算法。矩形求交也是最简单的求交算法之一,原理就是根据两个矩形的最大、最小坐标判断。图4展示了两个矩形没有交集的各种情况:
图4 矩形没有交集的几种情况
图5展示了两个矩形有交集的各种情况:
图5 矩形有交集的几种情况
从图4和图5可以分析出来两个矩形有交集的几何坐标规律,就是在x坐标方向和y坐标方向分别满足最大值最小值法则,简单解释这个法则就是每个矩形在 每个方向上的坐标最大值都要大于另一个矩形在这个坐标方向上的坐标最小值,否则在这个方向上就不能保证一定有位置重叠。由以上分析,判断两个矩形是否相交 的算法就可以如下实现:
186 bool IsRectIntersect(const Rect& rc1, const Rect& rc2) 187 { 188 return ( (std::max(rc1.p1.x, rc1.p2.x) >= std::min(rc2.p1.x, rc2.p2.x)) 189 && (std::max(rc2.p1.x, rc2.p2.x) >= std::min(rc1.p1.x, rc1.p2.x)) 190 && (std::max(rc1.p1.y, rc1.p2.y) >= std::min(rc2.p1.y, rc2.p2.y)) 191 && (std::max(rc2.p1.y, rc2.p2.y) >= std::min(rc1.p1.y, rc1.p2.y)) ); 192 }
完成了排斥试验和跨立试验的算法,最后判断直线段是否有交点的算法就水到渠成了:
204 bool IsLineSegmentIntersect(const LineSeg& ls1, const LineSeg& ls2) 205 { 206 if(IsLineSegmentExclusive(ls1, ls2)) //排斥实验 207 { 208 return false; 209 } 210 //( P1 - Q1 ) ×'a1?( Q2 - Q1 ) 211 double p1xq = CrossProduct(ls1.ps.x - ls2.ps.x, ls1.ps.y - ls2.ps.y, 212 ls2.pe.x - ls2.ps.x, ls2.pe.y - ls2.ps.y); 213 //( P2 - Q1 ) ×'a1?( Q2 - Q1 ) 214 double p2xq = CrossProduct(ls1.pe.x - ls2.ps.x, ls1.pe.y - ls2.ps.y, 215 ls2.pe.x - ls2.ps.x, ls2.pe.y - ls2.ps.y); 216 217 //( Q1 - P1 ) ×'a1?( P2 - P1 ) 218 double q1xp = CrossProduct(ls2.ps.x - ls1.ps.x, ls2.ps.y - ls1.ps.y, 219 ls1.pe.x - ls1.ps.x, ls1.pe.y - ls1.ps.y); 220 //( Q2 - P1 ) ×'a1?( P2 - P1 ) 221 double q2xp = CrossProduct(ls2.pe.x - ls1.ps.x, ls2.pe.y - ls1.ps.y, 222 ls1.pe.x - ls1.ps.x, ls1.pe.y - ls1.ps.y); 223 224 //跨立实验 225 return ( (p1xq * p2xq <= 0.0) && (q1xp * q2xp <= 0.0) ); 226 }
IsLineSegmentExclusive()函数就是调用IsRectIntersect()函数根据结果做排斥判断,此处不再列出代码。