GIS算法的計算幾何基礎(chǔ)

1、GIS算法的計算幾何基礎(chǔ)矢量的概念：如果一條線段的端點是有次序之分的，我們把這種線段成為有向線段(directed segment) 。如果有向線段p1p2 的起點 p1 在坐標原點，我們可以把它稱為矢量 (vector)p2 。矢量加減法：設二維矢量P = ( x1, y1 )， Q = ( x2 , y2 )，則矢量加法定義為：P + Q = ( x1 + x2 , y1 + y2 )，矢量減法定義為：P - Q = ( x1 - x2 , y1 - y2 )。顯然有性質(zhì)P + Q = Q + P ， P - Q = - ( Q - P )。矢量叉積 ：計算矢量叉積是與直線和線段相關(guān)算法的

2、核心部分。設矢量 P = ( x1, y1 )， Q = ( x2, y2 )，則矢量叉積定義為由(0,0)、pl、p2和p1+p2所組成的平行四邊形的帶符號的面積，即：P X Q = x1*y2 - x2*y1 ,其結(jié)果是一個標量。顯然有性質(zhì) P X Q = - ( Q X P )和 P X ( - Q ) = -( P X Q ) o兩點的加減法就是矢量相加減，而點的乘法則看作矢量叉積。叉積的一個非常重要性質(zhì)是可以通過它的符號判斷兩矢量相互之間的順逆時針關(guān)系：若P X Q > 0 ,則P在Q的順時針方向。若P X Q < 0 ,則P在Q的逆時針方向。若P X Q = 0 ,則P

3、與Q共線，但可能同向也可能反向。折線段的拐向判斷：折線段的拐向判斷方法可以直接由矢量叉積的性質(zhì)推出。對于有公共端點的線段p0p1和p1p2,通過計算(p2 - p0) X (p1 - p0)的符號便可以確定折線段的拐向：若(p2 - p0) X (p1 - p0) > 0,則p0p1在pl點拐向右側(cè)后得到p1p2o若(p2 - p0) X (p1 - p0) < 0,則p0p1在pl點拐向左側(cè)后得到p1p2o若(p2 - p0) X (p1-p0） = 0,則 p0、pl、p2 三點共線。具體情況可參照下圖:p2pO京，篙霍東能左g磊就I）？則tp2-plpO若如2 -例）X Cp

4、l -一 點共線.判斷點是否在線段上： 設點為Q線段為P1P2 ,判斷點Q在該線段上的依據(jù)是：（Q - P1 ） X （ P2 - P1 ） = 0 且Q在以P1 ,P2為對角頂點的矩形內(nèi)。前者保證Q點在直線P1P2上，后者是保證Q點不在線段P1P2的延長線或反向延長線上，對于這一步驟的判斷可以用以下過程實現(xiàn)：if min(Xp1,Xp2) <= Xq <= max(Xp1,Xp2) and min(Yp1,Yp2) <= Yq <= max(Yp1,Yp2) return true;elsereturn false;特別要注意的是，由于需要考慮水平線段和垂直線段兩種特

5、殊情況，min(xi,xj)<=xk<=max(xi,xj) 和 min(yi,yj)<=yk<=max(yi,yj) 兩個條件必須同 時滿足才能返回真值。判斷兩線段是否相交： 我們分兩步確定兩條線段是否相交：(1) 快速排斥試驗設以線段P1P2 為對角線的矩形為R， 設以線段Q1Q2 為對角線的矩形為T,如果R和T不相交，顯然兩線段不會相交。(2) 跨立試驗如果兩線段相交，則兩線段必然相互跨立對方。若P1P2跨立Q1Q2,則矢量(P1 - Q1 ) 和(P2 - Q1 ) 位于矢量(Q2 - Q1 ) 的兩側(cè)，即(P1 - Q1 ) X ( Q2 - Q1 ) * (

6、 P2 - Q1 ) X ( Q2 - Q1 ) < 0 。上式可改寫成(P1 - Q1 ) X ( Q2 - Q1 ) * ( Q2 - Q1 ) 乂 ( P2 - Q1 ) > 0 當(P1 - Q1 ) X ( Q2 - Q1 ) = 0 時，說明(P1 - Q1 ) 和(Q2 - Q1 ) 共線，但是因為已經(jīng)通過快速排斥試驗，所以P1 一定在線段Q1Q2上；同理，(Q2 - Q1 ) X (P2 - Q1 ) = 0 說明P2 一定在線段Q1Q2上。所以判斷 P1P2跨立 Q1Q2勺依據(jù)是：(P1 - Q1 ) X ( Q2 - Q1 ) * ( Q2 - Q1 ) X (

7、 P2 - Q1 ) >= 0 0同理判斷 Q1Q流立 P1P2的依據(jù)是：(Q1 - P1 ) X ( P2 - P1 ) * ( P2 - P1 ) X ( Q2 - P1 ) >= 0 。判斷線段和直線是否相交：如果線段P1P2和直線Q1Q才目交，貝U P1P2跨立Q1Q2即：(P1 - Q1 ) X ( Q2 - Q1 ) * ( Q2 - Q1 ) X ( P2 - Q1 ) >= 0 。判斷矩形是否包含點：只要判斷該點的橫坐標和縱坐標是否夾在矩形的左右邊和上下邊之間。判斷線段、折線、多邊形是否在矩形中：因為矩形是個凸集，所以只要判斷所有端點是否都在矩形中就可以了判斷

8、矩形是否在矩形中：只要比較左右邊界和上下邊界就可以了。判斷圓是否在矩形中：很容易證明，圓在矩形中的充要條件是：圓心在矩形中且圓的半徑小于等于圓心到矩形四邊的距離的最小值。判斷點是否在多邊形中：判斷點 P 是否在多邊形中是計算幾何中一個非常基本但是十分重要的算法。以點P為端點，向左方作射線L,由于多邊形是有界的，所以射線 L的左端一定在多邊形外，考慮沿著L 從無窮遠處開始自左向右移動，遇到和多邊形的第一個交點的時候，進入到了多邊形的內(nèi)部，遇到第二個交點的時候，離開了多邊形，所以很容易看出當L和多邊形的交點數(shù)目C是奇數(shù)的時候，P在多邊形內(nèi)，是偶數(shù)的話P 在多邊形外。但是有些特殊情況要加以考慮。如圖

9、下圖(a)(b)(c)(d) 所示。在圖(a)中， L 和多邊形的頂點相交，這時候交點只能計算一個；在圖(b) 中， L 和多邊形頂點的交點不應被計算；在圖(c) 和 (d) 中，L 和多邊形的一條邊重合，這條邊應該被忽略不計。如果L 和多邊形的一條邊重合，這條邊應該被忽略不計。為了統(tǒng)一起見，我們在計算射線 L和多邊形的交點的時候，1。對于多 邊形的水平邊不作考慮；2。對于多邊形的頂點和L相交的情況，如果該頂點是 其所屬的邊上縱坐標較大的頂點，則計數(shù)，否則忽略； 3。對于P在多邊形邊上 的情形，直接可判斷P屬于多邊行。由此得出算法的偽代碼如下：count - 0;以P為端點，作從右向左的射線L

10、;for多邊形的每條邊sdo if P 在邊s上then return true;if s 不是水平的then if s的一個端點在L上if 該端點是s 兩端點中縱坐標較大的端點then count count+1else if s 和 L 相交then count count+1;if count mod 2 = 1then return true;else return false;其中做射線L的方法是：設P'的縱坐標和P相同，橫坐標為正無窮大 （很大的一個正數(shù)），則P和P'就確定了射線L。判斷點是否在多邊形中的這個算法的時間復雜度為O（n）。另外還有一種算法是用帶符號的三

11、角形面積之和與多邊形面積進行比較，這種算法由于使用浮點數(shù)運算所以會帶來一定誤差，不推薦大家使用。判斷線段是否在多邊形內(nèi)：線段在多邊形內(nèi)的一個必要條件是線段的兩個端點都在多邊形內(nèi)，但由于多邊形可能為凹，所以這不能成為判斷的充分條件。如果線段和多邊形的某條邊內(nèi)交 （兩線段內(nèi)交是指兩線段相交且交點不在兩線段的端點）， 因為多邊形的邊的左右兩側(cè)分屬多邊形內(nèi)外不同部分，所以線段一定會有一部分在多邊形外（見圖a）。于是我們得到線段在多邊形內(nèi)的第二個必要條件：線段和多邊形的所有邊都不內(nèi)交。線段和多邊形交于線段的兩端點并不會影響線段是否在多邊形內(nèi)；但是如果多邊形的某個頂點和線段相交，還必須判斷兩相鄰交點之間的

12、線段是否包含于多邊形內(nèi)部（反例見圖b）。因此我們可以先求出所有和線段相交的多邊形的頂點，然后按照X-Y坐標排序（X坐標小的排在前面，對于X坐標相同的點，Y坐標小的排在前面，這 種排序準則也是為了保證水平和垂直情況的判斷正確），這樣相鄰的兩個點就是 在線段上相鄰的兩交點，如果任意相鄰兩點的中點也在多邊形內(nèi)，則該線段一定在多邊形內(nèi)。證明如下：命題1:如果線段和多邊形的兩相鄰交點 P1尸2的中點P'也在多邊形內(nèi)， 則P1, P2之間的所有點都在多邊形內(nèi)。證明：假設P1,P2之間含有不在多邊形內(nèi)的點，不妨設該點為Q,在P1, P' 之間，因為多邊形是閉合曲線，所以其內(nèi)外部之間有界，而P

13、1屬于多邊行內(nèi)部， Q屬于多邊性外部，P'屬于多邊性內(nèi)部，P1-Q-P'完全連續(xù)，所以P1CO QP'一定 跨越多邊形的邊界，因此在P1,P'之間至少還有兩個該線段和多邊形的交點，這 和P1P2是相鄰兩交點矛盾，故命題成立。證畢。由命題1直接可得出推論：推論2:設多邊形和線段PQ的交點依次為P1,P2,Pn,其中Pi和Pi+1 是相鄰兩交點，線段PQ在多邊形內(nèi)的充要條件是：P, Q在多邊形內(nèi)且對于i =1, 2,n -1 , Pi ,Pi+1 的中點也在多邊形內(nèi)。在實際編程中，沒有必要計算所有的交點，首先應判斷線段和多邊形的邊是否內(nèi)交，倘若線段和多邊形的某條邊內(nèi)

14、交則線段一定在多邊形外；如果線段和多邊形的每一條邊都不內(nèi)交，則線段和多邊形的交點一定是線段的端點或者多邊形的頂點，只要判斷點是否在線段上就可以了。至此我們得出算法如下：if線端PQ的端點不都在多邊形內(nèi)then return false;點集 pointSet 初始化為空;for 多邊形的每條邊sdo if 線段的某個端點在s 上then 將該端點加入pointSet;else if s的某個端點在線段PQ±then 將該端點加入pointSet;else if s 和線段PQ相交/這時候已經(jīng)可以肯定是內(nèi)交了then return false;將pointSet中的點按照X-Y坐標排序

15、;for pointSet 中每兩個相鄰點pointSeti , pointSet i+1do if pointSeti , pointSet i+1的中點不在多邊形中then return false;return true;這個過程中的排序因為交點數(shù)目肯定遠小于多邊形的頂點數(shù)目n， 所以最多是常數(shù)級的復雜度，幾乎可以忽略不計。因此算法的時間復雜度也是O(n)。判斷折線是否在多邊形內(nèi)：只要判斷折線的每條線段是否都在多邊形內(nèi)即可。設折線有 m條線段， 多邊形有n 個頂點，則該算法的時間復雜度為O(m*n)。判斷多邊形是否在多邊形內(nèi)：只要判斷多邊形的每條邊是否都在多邊形內(nèi)即可。判斷一個有m個頂點

16、的多邊形是否在一個有n 個頂點的多邊形內(nèi)復雜度為O(m*n)。判斷矩形是否在多邊形內(nèi)：將矩形轉(zhuǎn)化為多邊形，然后再判斷是否在多邊形內(nèi)。判斷圓是否在多邊形內(nèi)：只要計算圓心到多邊形的每條邊的最短距離，如果該距離大于等于圓半徑則該圓在多邊形內(nèi)。計算圓心到多邊形每條邊最短距離的算法在后文闡述。判斷點是否在圓內(nèi)：計算圓心到該點的距離，如果小于等于半徑則該點在圓內(nèi)。判斷線段、折線、矩形、多邊形是否在圓內(nèi)：因為圓是凸集，所以只要判斷是否每個頂點都在圓內(nèi)即可。判斷圓是否在圓內(nèi)：設兩圓為01,02半徑分別為r1, r2,要判斷02是否在01內(nèi)。先比較 ri , r2的大小，如果r1<r2則02不可能在01內(nèi)

17、；否則如果兩圓心的距離大于 r1 - r2 ，則02不在01 內(nèi)；否則02在01內(nèi)。計算點到線段的最近點：如果該線段平行于X軸(丫軸)，則過點point作該線段所在直線的垂線，垂足很容易求得，然后計算出垂足，如果垂足在線段上則返回垂足，否則返回離垂足近的端點；如果該線段不平行于X軸也不平行于丫軸，則斜率存在且不 為0。設線段的兩端點為pti和pt2 ,斜率為：k = ( pt2.y - pti. y ) / (pt2.x - pt1.x ); 該直線方程為：y = k* ( x - pt1.x) + pt1.y 。其垂線的斜率為- 1/ k ，垂線方程為：y = (-1/k) * (x - p

18、oint.x) + point.y。聯(lián)立兩直線方程解得：x = ( kA2 * pti.x + k * (point.y - pti.y )+ point.x ) / ( kA2 + 1), y = k * ( x - pti.x) + pti.y;然后再判斷垂足是否在線段上，如果在線段上則返回垂足；如果不在則計算兩端點到垂足的距離，選擇距離垂足較近的端點返回。計算點到折線、矩形、多邊形的最近點：只要分別計算點到每條線段的最近點，記錄最近距離，取其中最近距離最小的點即可。計算點到圓的最近距離及交點坐標如果該點在圓心，因為圓心到圓周任一點的距離相等，返回UNDEFINE。 D連接點P和圓心O,如

19、果PORT于X軸，則根據(jù)P在O的左邊還是右邊計算出最近點的橫坐標為centerPoint.x - radius 或 centerPoint.x +radius。如果PO平行于Y軸，則根據(jù)P在O的上邊還是下邊計算出最近點的縱坐標為 centerPoint.y -+radius 或 centerPoint.y - radius 。如果 PO不平行 于X軸和Y軸，則PO的斜率存在且不為0,這時直線POM率為k = ( P.y - O.y ) / ( P.x - O.x ) o直線PO的方程為：y = k * ( x - P.x) + P.y。設圓方程為：(x-O.x ) A2 + ( y - O.y

20、 ) A2 = r A2,聯(lián)立兩方程組可以解出直線 PO和圓的交點，取其中離P 點較近的交點即可。計算兩條共線的線段的交點：對于兩條共線的線段，它們之間的位置關(guān)系有下圖所示的幾種情況。圖(a) 中兩條線段沒有交點；圖 (b) 和 (d) 中兩條線段有無窮焦點；圖 (c) 中兩條線段有一個交點。設line1 是兩條線段中較長的一條，line2 是較短的一條，如果 line1 包含了 line2 的兩個端點，則是圖(d) 的情況，兩線段有無窮交點；如果 line1 只包含 line2 的一個端點，那么如果line1 的某個端點等于被line1包含的 line2 的那個端點，則是圖(c) 的情況，這

21、時兩線段只有一個交點，否則就是圖 (b) 的情況，兩線段也是有無窮的交點；如果line1 不包含 line2 的任何端點，則是圖(a) 的情況，這時兩線段沒有交點。計算線段或直線與線段的交點:設一條線段為L0 = P1P2,另一條線段或直線為L1 = Q1Q2,要計算的 就是 L0 和 L1 的交點。1首先判斷L0 和 L1 是否相交(方法已在前文討論過)，如果不相交則沒有交點，否則說明L0和L1 一定有交點，下面就將L0和L1都看作直線來考慮。2 .如果P1和P2橫坐標相同，即L0平行于Y軸a)若L1也平行于Y軸，i. 若 P1 的縱坐標和Q1 的縱坐標相同，說明L0 和 L1 共線，假如L

22、1 是直線的話他們有無窮的交點，假如L1 是線段的話可用" 計算兩條共線線段的交點 " 的算法求他們的交點(該方法在前文已討論過)；ii. 否則說明L0 和 L1 平行，他們沒有交點；b)若L1不平行于Y軸，則交點橫坐標為P1的橫坐標，代入到L1的直線方程中可以計算出交點縱坐標；3 .如果P1和P2橫坐標不同，但是Q1和Q2橫坐標相同，即L1平行于Y 軸， 則交點橫坐標為Q1 的橫坐標，代入到 L0 的直線方程中可以計算出交點縱坐標；4 .如果P1和P2縱坐標相同，即L0平行于X軸a）若L1也平行于X軸，i. 若 P1 的橫坐標和Q1 的橫坐標相同，說明L0 和 L1 共線

23、，假如L1 是直線的話他們有無窮的交點，假如L1 是線段的話可用" 計算兩條共線線段的交點 " 的算法求他們的交點（該方法在前文已討論過）；ii. 否則說明L0 和 L1 平行，他們沒有交點；b）若L1不平行于X軸，則交點縱坐標為P1的縱坐標，代入到L1的直 線方程中可以計算出交點橫坐標；5 .如果P1和P2縱坐標不同，但是Q1和Q2縱坐標相同，即L1平行于X 軸， 則交點縱坐標為Q1 的縱坐標，代入到 L0 的直線方程中可以計算出交點橫坐標；6 剩下的情況就是L1 和 L0 的斜率均存在且不為0 的情況a） 計算出 L0 的斜率K0， L1 的斜率 K1 ；b） 如果 K

24、1 = K2i . 如果 Q1 在 L0 上， 則說明 L0 和 L1 共線， 假如 L1 是直線的話有無窮交點，假如L1 是線段的話可用"計算兩條共線線段的交點"的算法求他們的交點（該方法在前文已討論過）；ii .如果Q1不在L0上，則說明L0和L1平行，他們沒有交點。c） 聯(lián)立兩直線的方程組可以解出交點來這個算法并不復雜，但是要分情況討論清楚，尤其是當兩條線段共線的情況需要單獨考慮，所以在前文將求兩條共線線段的算法單獨寫出來。另外， 一開始就先利用矢量叉乘判斷線段與線段（或直線）是否相交，如果結(jié)果是相交，那么在后面就可以將線段全部看作直線來考慮。需要注意的是，我們可以將

25、直線或線段方程改寫為ax+by+c=0的形式，這樣一來上述過程的部分步驟可以合并， 縮短了代碼長度，但是由于先要求出參數(shù)，這種算法將花費更多的時間。求線段或直線與折線、矩形、多邊形的交點：分別求與每條邊的交點即可。求線段或直線與圓的交點:設圓心為O,圓半徑為r,直線(或線段)L上的兩點為P1,P2。1 .如果L是線段且P1, P2都包含在圓O內(nèi)，則沒有交點；否則進行下 一步。2 .如果L平行于Y軸，a) 計算圓心到L 的距離 dis ；b) 如果 dis > r 則 L 和圓沒有交點；c) 利用勾股定理，可以求出兩交點坐標，但要注意考慮L 和圓的相切情況。3 .如果L平行于X軸，做法與L

26、平行于Y軸的情況類似；4 .如果L既不平行X軸也不平行Y軸，可以求出L的斜率K,然后列出 L 的點斜式方程，和圓方程聯(lián)立即可求解出L 和圓的兩個交點；5 . 如果 L 是線段，對于2， 3， 4 中求出的交點還要分別判斷是否屬于該線段的范圍內(nèi)。凸包的概念：點集Q的凸包(convex hull)是指一個最小凸多邊形，滿足 Q中的點或者在多邊形邊上或者在其內(nèi)。下圖中由紅色線段表示的多邊形就是點集Q=p0,p1,.p12 的凸包。凸包的求法：現(xiàn)在已經(jīng)證明了凸包算法的時間復雜度下界是 O(n*logn),但是當凸包 的頂點數(shù)h也被考慮進去的話，Krikpatrick 和Seidel的剪枝搜索算法可以達

27、 到O(n*logh),在漸進意義下達到最優(yōu)。最常用的凸包算法是 Graham掃描法和 Jarvis步進法。本文只簡單介紹一下 Graham掃描法，其正確性的證明和Jarvis 步進法的過程大家可以參考算法導論。對于一個有三個或以上點的點集 Q, Graham掃描法的過程如下：令p0為Q中Y-X坐標排序下最小的點S<p1,p2,.pm>為對其余點按以p0為中心的極角逆時針排序所得的 點集(如果有多個點有相同的極角，除了距p0最遠的點外全部移除壓p0進棧S壓p1進棧S壓p2進棧Sfor i 3 to mS的棧頂元素以及pi構(gòu)成的折線do while由S的棧頂元素的下一個元素、 段不拐

28、向左側(cè)對S彈棧壓pi進棧Sreturn S;此過程執(zhí)行后，棧S由底至頂?shù)脑鼐褪荙的凸包頂點按逆時針排列的 點序列。需要注意的是，我們對點按極角逆時針排序時，并不需要真正求出極角， 只需要求出任意兩點的次序就可以了。而這個步驟可以用前述的矢量叉積性質(zhì)實 現(xiàn)1984年圖靈獎獲得者，Pascal語言的發(fā)明人Niklaus Wirth教授曾經(jīng)給程序 下過一個經(jīng)典定義 程序=算法+數(shù)據(jù)結(jié)構(gòu)”，一針見血。在我看來，算法和數(shù)據(jù)結(jié) 構(gòu)就好象哲學中的方法論和認識論一樣，前者明確了解決一個問題的思路，而后者提供了解決問題的對象。算法”一詞來源于公元9世紀波斯數(shù)學家比阿勒.霍瓦里松的著作代數(shù)對話 錄，它是指完成一件事情的具體步驟和方法。因此說白了它并不神秘，計算機 中的算法，也不過就是將