用矢量的叉積判斷直線段是否有交
矢量叉積計算的另一個常用用途是直線段求交。求交算法是計算機圖形學的核心算法,也是體現速度和穩定性的重要標志,高效並且穩定的求交算法是任何一個 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()函數根據結果做排斥判斷,此處不再列出代碼。