給定頂點座標均是整點(或正方形格點)的簡單多邊形,皮克定理說明了其面積A和內(nèi)部格點數(shù)目i、邊上格點數(shù)目b的關(guān)系:A = i + b/2 - 1。
證明
因為所有簡單多邊形都可切割為一個三角形和另一個簡單多邊形。考慮一個簡單多邊形P,及跟P有一條共同邊的三角形T。若P符合皮克公式,則只要證明P加上T的PT亦符合皮克公式(I),與及三角形符合皮克公式(II),就可根據(jù)數(shù)學歸納法,對于所有簡單多邊形皮克公式都是成立的。
多邊形
設(shè)P和T的共同邊上有c個格點。
- P的面積: iP + bP/2 - 1
- T的面積: iT + bT/2 - 1
- PT的面積:
- (iT + iP + c - 2) + (bT- c + 2 + bP - c + 2 ) /2 - 1
- = iPT + bPT/2 - 1
三角形
證明分三部分:證明以下的圖形符合皮克定理:
- 所有平行于軸線的矩形;
- 以上述矩形的兩條鄰邊和對角線組成的直角三角形;
- 所有三角形(因為它們都可內(nèi)接于矩形內(nèi),將矩形分割成原三角形和至多3個第二點提到的直角三角形)。
矩形
設(shè)矩形R長邊短邊各有m,n個格點:
- AR = (m-1)(n-1)
- iR = (m-2)(n-2)
- bR = 2(m+n)-4
- iR + bR/2 - 1
- = (m-2)(n-2) + (m+n) - 2 - 1
- = mn - (m + n) +1
- = (m-1)(n-1)
直角三角形
易見兩條鄰邊和對角線組成的兩個直角三角形全等,且i,b相等。設(shè)其斜邊上有c個格點。
- b = m+n+c-3
- i = ((m-2)(n-2) - c + 2)/2
- i + b/2 - 1
- = ((m-2)(n-2) - c + 2)/2 + (m+n+c-3)/2 - 1
- = (m-2)(n-2)/2 + (m+n - 3)/2
- = (m-1)(n-1)/2
一般三角形
推廣
- 取格點的組成圖形的面積為一單位。在平行四邊形格點,皮克定理依然成立。套用于任意三角形格點,皮克定理則是A = 2i + b - 2。
- 對于非簡單的多邊形P,皮克定理A = i + b/2 - χ(P),其中χ(P)表示P的歐拉特征數(shù)。
- 高維推廣:Ehrhart多項式;一維:植樹問題。
- 皮克定理和歐拉公式(V-E+F=2)等價。
定理提出者
Georg Alexander Pick,1859年生于維也納,1943年死于特萊西恩施塔特集中營。
相關(guān)書籍
外部連結(jié)
en:Pick's theorem fr:Théorème de Pick it:Teorema di Pick pl:Wzór Picka ru:Теорема Пика
posted @
2009-10-06 17:59 wyiu 閱讀(144) |
評論 (0) |
編輯 收藏
學習了pick定理
給定頂點座標均是整點(或正方形格點)的簡單多邊形,皮克定理說明了其面積A和內(nèi)部格點數(shù)目i、邊上格點數(shù)目b的關(guān)系:A = i + b/2 - 1。
#include<iostream>
#include<math.h>
using namespace std;

#define MAXN 500


struct point
{int x,y;};
point p[MAXN];
int gcd(int a,int b)


{
if(b==0) return a;
else return gcd(b,a%b);
}
//計算多邊形面積,頂點按順時針或逆時針給出

double area_polygon(int n,point* p)
{
double s1=0,s2=0;
int i;
for (i=0;i<n;i++)
s1+=p[(i+1)%n].y*p[i].x,s2+=p[(i+1)%n].y*p[(i+2)%n].x;
return fabs(s1-s2)/2;
}

int main()


{
int tc,m,I,E,dx,dy,i,j,k;
scanf("%d",&tc);
for(k=1;k<=tc;k++)

{
scanf("%d",&m);
p[0].x=0;
p[0].y=0;
E=0;
for(i=1;i<=m;i++)

{
scanf("%d%d",&dx,&dy);
p[i].x=p[i-1].x+dx;
p[i].y=p[i-1].y+dy;
E+=gcd(abs(dx),abs(dy));
}
double area=area_polygon(m+1,p);
I=(2*area+2-E)/2;
printf("Scenario #%d:\n%d %d %.1lf\n\n",k,I,E,area);
}
return 0;
}

posted @
2009-10-06 17:58 wyiu 閱讀(241) |
評論 (0) |
編輯 收藏
一、引言
計算機的出現(xiàn)使得很多原本十分繁瑣的工作得以大幅度簡化,但是也有一些在人們直觀看來很容易的問題卻需要拿出一套并不簡單的通用解決方案,比如幾何問題。作為計算機科學的一個分支,計算幾何主要研究解決幾何問題的算法。在現(xiàn)代工程和數(shù)學領(lǐng)域,計算幾何在圖形學、機器人技術(shù)、超大規(guī)模集成電路設(shè)計和統(tǒng)計等諸多領(lǐng)域有著十分重要的應(yīng)用。在本文中,我們將對計算幾何常用的基本算法做一個全面的介紹,希望對您了解并應(yīng)用計算幾何的知識解決問題起到幫助。
二、目錄 本文整理的計算幾何基本概念和常用算法包括如下內(nèi)容:
矢量的概念
矢量加減法
矢量叉積
折線段的拐向判斷
判斷點是否在線段上
判斷兩線段是否相交
判斷線段和直線是否相交
判斷矩形是否包含點
判斷線段、折線、多邊形是否在矩形中
判斷矩形是否在矩形中
判斷圓是否在矩形中
判斷點是否在多邊形中
判斷線段是否在多邊形內(nèi)
判斷折線是否在多邊形內(nèi)
判斷多邊形是否在多邊形內(nèi)
判斷矩形是否在多邊形內(nèi)
判斷圓是否在多邊形內(nèi)
判斷點是否在圓內(nèi)
判斷線段、折線、矩形、多邊形是否在圓內(nèi)
判斷圓是否在圓內(nèi)
計算點到線段的最近點
計算點到折線、矩形、多邊形的最近點
計算點到圓的最近距離及交點坐標
計算兩條共線的線段的交點
計算線段或直線與線段的交點
求線段或直線與折線、矩形、多邊形的交點
求線段或直線與圓的交點
凸包的概念
凸包的求法
三、算法介紹矢量的概念:如果一條線段的端點是有次序之分的,我們把這種線段成為有向線段(directed segment)。如果有向線段p1p2的起點p1在坐標原點,我們可以把它稱為矢量(vector)p2。
矢量叉積:計算矢量叉積是與直線和線段相關(guān)算法的核心部分。設(shè)矢量P = (x1,y1) ,Q = (x2,y2),則矢量叉積定義為由(0,0)、p1、p2和p1+p2所組成的平行四邊形的帶符號的面積,即:P × Q = x1*y2 - x2*y1,其結(jié)果是一個標量。顯然有性質(zhì) P × Q = - ( Q × P ) 和 P × ( - Q ) = - ( P × Q )。一般在不加說明的情況下,本文下述算法中所有的點都看作矢量,兩點的加減法就是矢量相加減,而點的乘法則看作矢量叉積。
叉積的一個非常重要性質(zhì)是可以通過它的符號判斷兩矢量相互之間的順逆時針關(guān)系:
若 P × Q > 0 , 則P在Q的順時針方向。
若 P × Q < 0 , 則P在Q的逆時針方向。
若 P × Q = 0 , 則P與Q共線,但可能同向也可能反向。
折線段的拐向判斷: 折線段的拐向判斷方法可以直接由矢量叉積的性質(zhì)推出。對于有公共端點的線段p0p1和p1p2,通過計算(p2 - p0) × (p1 - p0)的符號便可以確定折線段的拐向:
若(p2 - p0) × (p1 - p0) > 0,則p0p1在p1點拐向右側(cè)后得到p1p2。
若(p2 - p0) × (p1 - p0) < 0,則p0p1在p1點拐向左側(cè)后得到p1p2。
若(p2 - p0) × (p1 - p0) = 0,則p0、p1、p2三點共線。
具體情況可參照下圖:

判斷點是否在線段上:
設(shè)點為Q,線段為P1P2 ,判斷點Q在該線段上的依據(jù)是:( Q - P1 ) × ( P2 - P1 ) = 0 且 Q 在以 P1,P2為對角頂點的矩形內(nèi)。前者保證Q點在直線P1P2上,后者是保證Q點不在線段P1P2的延長線或反向延長線上,對于這一步驟的判斷可以用以下過程實現(xiàn):
ON-SEGMENT(pi,pj,pk)
if min(xi,xj)<=xk<=max(xi,xj) and min(yi,yj)<=yk<=max(yi,yj)
then return true;
else return false;
特別要注意的是,由于需要考慮水平線段和垂直線段兩種特殊情況,min(xi,xj)<=xk<=max(xi,xj)和min(yi,yj)<=yk<=max(yi,yj)兩個條件必須同時滿足才能返回真值。
判斷兩線段是否相交:
我們分兩步確定兩條線段是否相交:
(1)快速排斥試驗
設(shè)以線段 P1P2 為對角線的矩形為R, 設(shè)以線段 Q1Q2 為對角線的矩形為T,如果R和T不相交,顯然兩線段不會相交。
(2)跨立試驗
如果兩線段相交,則兩線段必然相互跨立對方。若P1P2跨立Q1Q2 ,則矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的兩側(cè),即( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0。上式可改寫成( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0。當 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 時,說明 ( P1 - Q1 ) 和 ( Q2 - Q1 )共線,但是因為已經(jīng)通過快速排斥試驗,所以 P1 一定在線段 Q1Q2上;同理,( Q2 - Q1 ) ×(P2 - Q1 ) = 0 說明 P2 一定在線段 Q1Q2上。所以判斷P1P2跨立Q1Q2的依據(jù)是:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。同理判斷Q1Q2跨立P1P2的依據(jù)是:( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。具體情況如下圖所示:

在相同的原理下,對此算法的具體的實現(xiàn)細節(jié)可能會與此有所不同,除了這種過程外,大家也可以參考《算法導(dǎo)論》上的實現(xiàn)。
判斷線段和直線是否相交:
有了上面的基礎(chǔ),這個算法就很容易了。如果線段P1P2和直線Q1Q2相交,則P1P2跨立Q1Q2,即:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
判斷矩形是否包含點:
只要判斷該點的橫坐標和縱坐標是否夾在矩形的左右邊和上下邊之間。
判斷線段、折線、多邊形是否在矩形中:
因為矩形是個凸集,所以只要判斷所有端點是否都在矩形中就可以了。
判斷矩形是否在矩形中:
只要比較左右邊界和上下邊界就可以了。
判斷圓是否在矩形中:
很容易證明,圓在矩形中的充要條件是:圓心在矩形中且圓的半徑小于等于圓心到矩形四邊的距離的最小值。
判斷點是否在多邊形中:
判斷點P是否在多邊形中是計算幾何中一個非常基本但是十分重要的算法。以點P為端點,向左方作射線L,由于多邊形是有界的,所以射線L的左端一定在多邊形外,考慮沿著L從無窮遠處開始自左向右移動,遇到和多邊形的第一個交點的時候,進入到了多邊形的內(nèi)部,遇到第二個交點的時候,離開了多邊形,……所以很容易看出當L和多邊形的交點數(shù)目C是奇數(shù)的時候,P在多邊形內(nèi),是偶數(shù)的話P在多邊形外。
但是有些特殊情況要加以考慮。如圖下圖(a)(b)(c)(d)所示。在圖(a)中,L和多邊形的頂點相交,這時候交點只能計算一個;在圖(b)中,L和多邊形頂點的交點不應(yīng)被計算;在圖(c)和(d) 中,L和多邊形的一條邊重合,這條邊應(yīng)該被忽略不計。如果L和多邊形的一條邊重合,這條邊應(yīng)該被忽略不計。

為了統(tǒng)一起見,我們在計算射線L和多邊形的交點的時候,1。對于多邊形的水平邊不作考慮;2。對于多邊形的頂點和L相交的情況,如果該頂點是其所屬的邊上縱坐標較大的頂點,則計數(shù),否則忽略;3。對于P在多邊形邊上的情形,直接可判斷P屬于多邊行。由此得出算法的偽代碼如下:
count ← 0;
以P為端點,作從右向左的射線L;
for 多邊形的每條邊s
do if P在邊s上
then return true;
if s不是水平的
then if s的一個端點在L上
if 該端點是s兩端點中縱坐標較大的端點
then count ← count+1
else if s和L相交
then count ← count+1;
if count mod 2 = 1
then return true;
else return false;
其中做射線L的方法是:設(shè)P'的縱坐標和P相同,橫坐標為正無窮大(很大的一個正數(shù)),則P和P'就確定了射線L。
判斷點是否在多邊形中的這個算法的時間復(fù)雜度為O(n)。
另外還有一種算法是用帶符號的三角形面積之和與多邊形面積進行比較,這種算法由于使用浮點數(shù)運算所以會帶來一定誤差,不推薦大家使用。
判斷線段是否在多邊形內(nèi):
線段在多邊形內(nèi)的一個必要條件是線段的兩個端點都在多邊形內(nèi),但由于多邊形可能為凹,所以這不能成為判斷的充分條件。如果線段和多邊形的某條邊內(nèi)交(兩線段內(nèi)交是指兩線段相交且交點不在兩線段的端點),因為多邊形的邊的左右兩側(cè)分屬多邊形內(nèi)外不同部分,所以線段一定會有一部分在多邊形外(見圖a)。于是我們得到線段在多邊形內(nèi)的第二個必要條件:線段和多邊形的所有邊都不內(nèi)交。
線段和多邊形交于線段的兩端點并不會影響線段是否在多邊形內(nèi);但是如果多邊形的某個頂點和線段相交,還必須判斷兩相鄰交點之間的線段是否包含于多邊形內(nèi)部(反例見圖b)。

因此我們可以先求出所有和線段相交的多邊形的頂點,然后按照X-Y坐標排序(X坐標小的排在前面,對于X坐標相同的點,Y坐標小的排在前面,這種排序準則也是為了保證水平和垂直情況的判斷正確),這樣相鄰的兩個點就是在線段上相鄰的兩交點,如果任意相鄰兩點的中點也在多邊形內(nèi),則該線段一定在多邊形內(nèi)。
證明如下:
命題1:
如果線段和多邊形的兩相鄰交點P1 ,P2的中點P' 也在多邊形內(nèi),則P1, P2之間的所有點都在多邊形內(nèi)。
證明:
假設(shè)P1,P2之間含有不在多邊形內(nèi)的點,不妨設(shè)該點為Q,在P1, P'之間,因為多邊形是閉合曲線,所以其內(nèi)外部之間有界,而P1屬于多邊行內(nèi)部,Q屬于多邊性外部,P'屬于多邊性內(nèi)部,P1-Q-P'完全連續(xù),所以P1Q和QP'一定跨越多邊形的邊界,因此在P1,P'之間至少還有兩個該線段和多邊形的交點,這和P1P2是相鄰兩交點矛盾,故命題成立。證畢。
由命題1直接可得出推論:
推論2:
設(shè)多邊形和線段PQ的交點依次為P1,P2,……Pn,其中Pi和Pi+1是相鄰兩交點,線段PQ在多邊形內(nèi)的充要條件是:P,Q在多邊形內(nèi)且對于i =1, 2,……, n-1,Pi ,Pi+1的中點也在多邊形內(nèi)。
在實際編程中,沒有必要計算所有的交點,首先應(yīng)判斷線段和多邊形的邊是否內(nèi)交,倘若線段和多邊形的某條邊內(nèi)交則線段一定在多邊形外;如果線段和多邊形的每一條邊都不內(nèi)交,則線段和多邊形的交點一定是線段的端點或者多邊形的頂點,只要判斷點是否在線段上就可以了。
至此我們得出算法如下:
if 線端PQ的端點不都在多邊形內(nèi)
then return false;
點集pointSet初始化為空;
for 多邊形的每條邊s
do if 線段的某個端點在s上
then 將該端點加入pointSet;
else if s的某個端點在線段PQ上
then 將該端點加入pointSet;
else if s和線段PQ相交 // 這時候已經(jīng)可以肯定是內(nèi)交了
then return false;
將pointSet中的點按照X-Y坐標排序;
for pointSet中每兩個相鄰點 pointSet[i] , pointSet[ i+1]
do if pointSet[i] , pointSet[ i+1] 的中點不在多邊形中
then return false;
return true;
這個過程中的排序因為交點數(shù)目肯定遠小于多邊形的頂點數(shù)目n,所以最多是常數(shù)級的復(fù)雜度,幾乎可以忽略不計。因此算法的時間復(fù)雜度也是O(n)。
判斷折線是否在多邊形內(nèi):
只要判斷折線的每條線段是否都在多邊形內(nèi)即可。設(shè)折線有m條線段,多邊形有n個頂點,則該算法的時間復(fù)雜度為O(m*n)。
判斷多邊形是否在多邊形內(nèi):
只要判斷多邊形的每條邊是否都在多邊形內(nèi)即可。判斷一個有m個頂點的多邊形是否在一個有n個頂點的多邊形內(nèi)復(fù)雜度為O(m*n)。
判斷矩形是否在多邊形內(nèi):
將矩形轉(zhuǎn)化為多邊形,然后再判斷是否在多邊形內(nèi)。
判斷圓是否在多邊形內(nèi):
只要計算圓心到多邊形的每條邊的最短距離,如果該距離大于等于圓半徑則該圓在多邊形內(nèi)。計算圓心到多邊形每條邊最短距離的算法在后文闡述。
判斷點是否在圓內(nèi):
計算圓心到該點的距離,如果小于等于半徑則該點在圓內(nèi)。
判斷線段、折線、矩形、多邊形是否在圓內(nèi):
因為圓是凸集,所以只要判斷是否每個頂點都在圓內(nèi)即可。
判斷圓是否在圓內(nèi):
設(shè)兩圓為O1,O2,半徑分別為r1, r2,要判斷O2是否在O1內(nèi)。先比較r1,r2的大小,如果r1<r2則O2不可能在O1內(nèi);否則如果兩圓心的距離大于r1 - r2 ,則O2不在O1內(nèi);否則O2在O1內(nèi)。
計算點到線段的最近點:
如果該線段平行于X軸(Y軸),則過點point作該線段所在直線的垂線,垂足很容易求得,然后計算出垂足,如果垂足在線段上則返回垂足,否則返回離垂足近的端點;如果該線段不平行于X軸也不平行于Y軸,則斜率存在且不為0。設(shè)線段的兩端點為pt1和pt2,斜率為:k = ( pt2.y - pt1. y ) / (pt2.x - pt1.x );該直線方程為:y = k* ( x - pt1.x) + pt1.y。其垂線的斜率為 - 1 / k,垂線方程為:y = (-1/k) * (x - point.x) + point.y 。
聯(lián)立兩直線方程解得:x = ( k^2 * pt1.x + k * (point.y - pt1.y ) + point.x ) / ( k^2 + 1) ,y = k * ( x - pt1.x) + pt1.y;然后再判斷垂足是否在線段上,如果在線段上則返回垂足;如果不在則計算兩端點到垂足的距離,選擇距離垂足較近的端點返回。
計算點到折線、矩形、多邊形的最近點:
只要分別計算點到每條線段的最近點,記錄最近距離,取其中最近距離最小的點即可。
計算點到圓的最近距離及交點坐標:
如果該點在圓心,因為圓心到圓周任一點的距離相等,返回UNDEFINED。
連接點P和圓心O,如果PO平行于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,這時直線PO斜率為k = ( P.y - O.y )/ ( P.x - O.x )。直線PO的方程為:y = k * ( x - P.x) + P.y。設(shè)圓方程為:(x - O.x ) ^2 + ( y - O.y ) ^2 = r ^2,聯(lián)立兩方程組可以解出直線PO和圓的交點,取其中離P點較近的交點即可。
計算兩條共線的線段的交點:
對于兩條共線的線段,它們之間的位置關(guān)系有下圖所示的幾種情況。圖(a)中兩條線段沒有交點;圖 (b) 和 (d) 中兩條線段有無窮焦點;圖 (c) 中兩條線段有一個交點。設(shè)line1是兩條線段中較長的一條,line2是較短的一條,如果line1包含了line2的兩個端點,則是圖(d)的情況,兩線段有無窮交點;如果line1只包含line2的一個端點,那么如果line1的某個端點等于被line1包含的line2的那個端點,則是圖(c)的情況,這時兩線段只有一個交點,否則就是圖(b)的情況,兩線段也是有無窮的交點;如果line1不包含line2的任何端點,則是圖(a)的情況,這時兩線段沒有交點。

計算線段或直線與線段的交點:
設(shè)一條線段為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共線,假如L1是直線的話他們有無窮的交點,假如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共線,假如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) 如果K1 = K2
i. 如果Q1在L0上,則說明L0和L1共線,假如L1是直線的話有無窮交點,假如L1是線段的話可用"計算兩條共線線段的交點"的算法求他們的交點(該方法在前文已討論過);
ii. 如果Q1不在L0上,則說明L0和L1平行,他們沒有交點。
c) 聯(lián)立兩直線的方程組可以解出交點來
這個算法并不復(fù)雜,但是要分情況討論清楚,尤其是當兩條線段共線的情況需要單獨考慮,所以在前文將求兩條共線線段的算法單獨寫出來。另外,一開始就先利用矢量叉乘判斷線段與線段(或直線)是否相交,如果結(jié)果是相交,那么在后面就可以將線段全部看作直線來考慮。需要注意的是,我們可以將直線或線段方程改寫為ax+by+c=0的形式,這樣一來上述過程的部分步驟可以合并,縮短了代碼長度,但是由于先要求出參數(shù),這種算法將花費更多的時間。
求線段或直線與折線、矩形、多邊形的交點:
分別求與每條邊的交點即可。
求線段或直線與圓的交點:
設(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平行于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)證明了凸包算法的時間復(fù)雜度下界是O(n*logn),但是當凸包的頂點數(shù)h也被考慮進去的話,Krikpatrick和Seidel的剪枝搜索算法可以達到O(n*logh),在漸進意義下達到最優(yōu)。最常用的凸包算法是Graham掃描法和Jarvis步進法。本文只簡單介紹一下Graham掃描法,其正確性的證明和Jarvis步進法的過程大家可以參考《算法導(dǎo)論》。
對于一個有三個或以上點的點集Q,Graham掃描法的過程如下:
令p0為Q中Y-X坐標排序下最小的點
設(shè)<p1,p2,...pm>為對其余點按以p0為中心的極角逆時針排序所得的點集(如果有多個點有相同的極角,除了距p0最遠的點外全部移除
壓p0進棧S
壓p1進棧S
壓p2進棧S
for i ← 3 to m
do while 由S的棧頂元素的下一個元素、S的棧頂元素以及pi構(gòu)成的折線段不拐向左側(cè)
對S彈棧
壓pi進棧S
return S;
此過程執(zhí)行后,棧S由底至頂?shù)脑鼐褪荙的凸包頂點按逆時針排列的點序列。需要注意的是,我們對點按極角逆時針排序時,并不需要真正求出極角,只需要求出任意兩點的次序就可以了。而這個步驟可以用前述的矢量叉積性質(zhì)實現(xiàn)。
四、結(jié)語
盡管人類對幾何學的研究從古代起便沒有中斷過,但是具體到借助計算機來解決幾何問題的研究,還只是停留在一個初級階段,無論從應(yīng)用領(lǐng)域還是發(fā)展前景來看,計算幾何學都值得我們認真學習、加以運用,希望這篇文章能帶你走進這個豐富多彩的世界。
posted @
2009-10-06 14:38 wyiu 閱讀(156) |
評論 (0) |
編輯 收藏
一、引言
計算機的出現(xiàn)使得很多原本十分繁瑣的工作得以大幅度簡化,但是也有一些在人們直觀看來很容易的問題卻需要拿出一套并不簡單的通用解決方案,比如幾何問題。作為計算機科學的一個分支,計算幾何主要研究解決幾何問題的算法。在現(xiàn)代工程和數(shù)學領(lǐng)域,計算幾何在圖形學、機器人技術(shù)、超大規(guī)模集成電路設(shè)計和統(tǒng)計等諸多領(lǐng)域有著十分重要的應(yīng)用。在本文中,我們將對計算幾何常用的基本算法做一個全面的介紹,希望對您了解并應(yīng)用計算幾何的知識解決問題起到幫助。
二、目錄 本文整理的計算幾何基本概念和常用算法包括如下內(nèi)容:
矢量的概念
矢量加減法
矢量叉積
折線段的拐向判斷
判斷點是否在線段上
判斷兩線段是否相交
判斷線段和直線是否相交
判斷矩形是否包含點
判斷線段、折線、多邊形是否在矩形中
判斷矩形是否在矩形中
判斷圓是否在矩形中
判斷點是否在多邊形中
判斷線段是否在多邊形內(nèi)
判斷折線是否在多邊形內(nèi)
判斷多邊形是否在多邊形內(nèi)
判斷矩形是否在多邊形內(nèi)
判斷圓是否在多邊形內(nèi)
判斷點是否在圓內(nèi)
判斷線段、折線、矩形、多邊形是否在圓內(nèi)
判斷圓是否在圓內(nèi)
計算點到線段的最近點
計算點到折線、矩形、多邊形的最近點
計算點到圓的最近距離及交點坐標
計算兩條共線的線段的交點
計算線段或直線與線段的交點
求線段或直線與折線、矩形、多邊形的交點
求線段或直線與圓的交點
凸包的概念
凸包的求法
三、算法介紹矢量的概念:如果一條線段的端點是有次序之分的,我們把這種線段成為有向線段(directed segment)。如果有向線段p1p2的起點p1在坐標原點,我們可以把它稱為矢量(vector)p2。
矢量叉積:計算矢量叉積是與直線和線段相關(guān)算法的核心部分。設(shè)矢量P = (x1,y1) ,Q = (x2,y2),則矢量叉積定義為由(0,0)、p1、p2和p1+p2所組成的平行四邊形的帶符號的面積,即:P × Q = x1*y2 - x2*y1,其結(jié)果是一個標量。顯然有性質(zhì) P × Q = - ( Q × P ) 和 P × ( - Q ) = - ( P × Q )。一般在不加說明的情況下,本文下述算法中所有的點都看作矢量,兩點的加減法就是矢量相加減,而點的乘法則看作矢量叉積。
叉積的一個非常重要性質(zhì)是可以通過它的符號判斷兩矢量相互之間的順逆時針關(guān)系:
若 P × Q > 0 , 則P在Q的順時針方向。
若 P × Q < 0 , 則P在Q的逆時針方向。
若 P × Q = 0 , 則P與Q共線,但可能同向也可能反向。
折線段的拐向判斷: 折線段的拐向判斷方法可以直接由矢量叉積的性質(zhì)推出。對于有公共端點的線段p0p1和p1p2,通過計算(p2 - p0) × (p1 - p0)的符號便可以確定折線段的拐向:
若(p2 - p0) × (p1 - p0) > 0,則p0p1在p1點拐向右側(cè)后得到p1p2。
若(p2 - p0) × (p1 - p0) < 0,則p0p1在p1點拐向左側(cè)后得到p1p2。
若(p2 - p0) × (p1 - p0) = 0,則p0、p1、p2三點共線。
具體情況可參照下圖:

判斷點是否在線段上:
設(shè)點為Q,線段為P1P2 ,判斷點Q在該線段上的依據(jù)是:( Q - P1 ) × ( P2 - P1 ) = 0 且 Q 在以 P1,P2為對角頂點的矩形內(nèi)。前者保證Q點在直線P1P2上,后者是保證Q點不在線段P1P2的延長線或反向延長線上,對于這一步驟的判斷可以用以下過程實現(xiàn):
ON-SEGMENT(pi,pj,pk)
if min(xi,xj)<=xk<=max(xi,xj) and min(yi,yj)<=yk<=max(yi,yj)
then return true;
else return false;
特別要注意的是,由于需要考慮水平線段和垂直線段兩種特殊情況,min(xi,xj)<=xk<=max(xi,xj)和min(yi,yj)<=yk<=max(yi,yj)兩個條件必須同時滿足才能返回真值。
判斷兩線段是否相交:
我們分兩步確定兩條線段是否相交:
(1)快速排斥試驗
設(shè)以線段 P1P2 為對角線的矩形為R, 設(shè)以線段 Q1Q2 為對角線的矩形為T,如果R和T不相交,顯然兩線段不會相交。
(2)跨立試驗
如果兩線段相交,則兩線段必然相互跨立對方。若P1P2跨立Q1Q2 ,則矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的兩側(cè),即( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0。上式可改寫成( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0。當 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 時,說明 ( P1 - Q1 ) 和 ( Q2 - Q1 )共線,但是因為已經(jīng)通過快速排斥試驗,所以 P1 一定在線段 Q1Q2上;同理,( Q2 - Q1 ) ×(P2 - Q1 ) = 0 說明 P2 一定在線段 Q1Q2上。所以判斷P1P2跨立Q1Q2的依據(jù)是:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。同理判斷Q1Q2跨立P1P2的依據(jù)是:( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。具體情況如下圖所示:

在相同的原理下,對此算法的具體的實現(xiàn)細節(jié)可能會與此有所不同,除了這種過程外,大家也可以參考《算法導(dǎo)論》上的實現(xiàn)。
判斷線段和直線是否相交:
有了上面的基礎(chǔ),這個算法就很容易了。如果線段P1P2和直線Q1Q2相交,則P1P2跨立Q1Q2,即:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
判斷矩形是否包含點:
只要判斷該點的橫坐標和縱坐標是否夾在矩形的左右邊和上下邊之間。
判斷線段、折線、多邊形是否在矩形中:
因為矩形是個凸集,所以只要判斷所有端點是否都在矩形中就可以了。
判斷矩形是否在矩形中:
只要比較左右邊界和上下邊界就可以了。
判斷圓是否在矩形中:
很容易證明,圓在矩形中的充要條件是:圓心在矩形中且圓的半徑小于等于圓心到矩形四邊的距離的最小值。
判斷點是否在多邊形中:
判斷點P是否在多邊形中是計算幾何中一個非常基本但是十分重要的算法。以點P為端點,向左方作射線L,由于多邊形是有界的,所以射線L的左端一定在多邊形外,考慮沿著L從無窮遠處開始自左向右移動,遇到和多邊形的第一個交點的時候,進入到了多邊形的內(nèi)部,遇到第二個交點的時候,離開了多邊形,……所以很容易看出當L和多邊形的交點數(shù)目C是奇數(shù)的時候,P在多邊形內(nèi),是偶數(shù)的話P在多邊形外。
但是有些特殊情況要加以考慮。如圖下圖(a)(b)(c)(d)所示。在圖(a)中,L和多邊形的頂點相交,這時候交點只能計算一個;在圖(b)中,L和多邊形頂點的交點不應(yīng)被計算;在圖(c)和(d) 中,L和多邊形的一條邊重合,這條邊應(yīng)該被忽略不計。如果L和多邊形的一條邊重合,這條邊應(yīng)該被忽略不計。

為了統(tǒng)一起見,我們在計算射線L和多邊形的交點的時候,1。對于多邊形的水平邊不作考慮;2。對于多邊形的頂點和L相交的情況,如果該頂點是其所屬的邊上縱坐標較大的頂點,則計數(shù),否則忽略;3。對于P在多邊形邊上的情形,直接可判斷P屬于多邊行。由此得出算法的偽代碼如下:
count ← 0;
以P為端點,作從右向左的射線L;
for 多邊形的每條邊s
do if P在邊s上
then return true;
if s不是水平的
then if s的一個端點在L上
if 該端點是s兩端點中縱坐標較大的端點
then count ← count+1
else if s和L相交
then count ← count+1;
if count mod 2 = 1
then return true;
else return false;
其中做射線L的方法是:設(shè)P'的縱坐標和P相同,橫坐標為正無窮大(很大的一個正數(shù)),則P和P'就確定了射線L。
判斷點是否在多邊形中的這個算法的時間復(fù)雜度為O(n)。
另外還有一種算法是用帶符號的三角形面積之和與多邊形面積進行比較,這種算法由于使用浮點數(shù)運算所以會帶來一定誤差,不推薦大家使用。
判斷線段是否在多邊形內(nèi):
線段在多邊形內(nèi)的一個必要條件是線段的兩個端點都在多邊形內(nèi),但由于多邊形可能為凹,所以這不能成為判斷的充分條件。如果線段和多邊形的某條邊內(nèi)交(兩線段內(nèi)交是指兩線段相交且交點不在兩線段的端點),因為多邊形的邊的左右兩側(cè)分屬多邊形內(nèi)外不同部分,所以線段一定會有一部分在多邊形外(見圖a)。于是我們得到線段在多邊形內(nèi)的第二個必要條件:線段和多邊形的所有邊都不內(nèi)交。
線段和多邊形交于線段的兩端點并不會影響線段是否在多邊形內(nèi);但是如果多邊形的某個頂點和線段相交,還必須判斷兩相鄰交點之間的線段是否包含于多邊形內(nèi)部(反例見圖b)。

因此我們可以先求出所有和線段相交的多邊形的頂點,然后按照X-Y坐標排序(X坐標小的排在前面,對于X坐標相同的點,Y坐標小的排在前面,這種排序準則也是為了保證水平和垂直情況的判斷正確),這樣相鄰的兩個點就是在線段上相鄰的兩交點,如果任意相鄰兩點的中點也在多邊形內(nèi),則該線段一定在多邊形內(nèi)。
證明如下:
命題1:
如果線段和多邊形的兩相鄰交點P1 ,P2的中點P' 也在多邊形內(nèi),則P1, P2之間的所有點都在多邊形內(nèi)。
證明:
假設(shè)P1,P2之間含有不在多邊形內(nèi)的點,不妨設(shè)該點為Q,在P1, P'之間,因為多邊形是閉合曲線,所以其內(nèi)外部之間有界,而P1屬于多邊行內(nèi)部,Q屬于多邊性外部,P'屬于多邊性內(nèi)部,P1-Q-P'完全連續(xù),所以P1Q和QP'一定跨越多邊形的邊界,因此在P1,P'之間至少還有兩個該線段和多邊形的交點,這和P1P2是相鄰兩交點矛盾,故命題成立。證畢。
由命題1直接可得出推論:
推論2:
設(shè)多邊形和線段PQ的交點依次為P1,P2,……Pn,其中Pi和Pi+1是相鄰兩交點,線段PQ在多邊形內(nèi)的充要條件是:P,Q在多邊形內(nèi)且對于i =1, 2,……, n-1,Pi ,Pi+1的中點也在多邊形內(nèi)。
在實際編程中,沒有必要計算所有的交點,首先應(yīng)判斷線段和多邊形的邊是否內(nèi)交,倘若線段和多邊形的某條邊內(nèi)交則線段一定在多邊形外;如果線段和多邊形的每一條邊都不內(nèi)交,則線段和多邊形的交點一定是線段的端點或者多邊形的頂點,只要判斷點是否在線段上就可以了。
至此我們得出算法如下:
if 線端PQ的端點不都在多邊形內(nèi)
then return false;
點集pointSet初始化為空;
for 多邊形的每條邊s
do if 線段的某個端點在s上
then 將該端點加入pointSet;
else if s的某個端點在線段PQ上
then 將該端點加入pointSet;
else if s和線段PQ相交 // 這時候已經(jīng)可以肯定是內(nèi)交了
then return false;
將pointSet中的點按照X-Y坐標排序;
for pointSet中每兩個相鄰點 pointSet[i] , pointSet[ i+1]
do if pointSet[i] , pointSet[ i+1] 的中點不在多邊形中
then return false;
return true;
這個過程中的排序因為交點數(shù)目肯定遠小于多邊形的頂點數(shù)目n,所以最多是常數(shù)級的復(fù)雜度,幾乎可以忽略不計。因此算法的時間復(fù)雜度也是O(n)。
判斷折線是否在多邊形內(nèi):
只要判斷折線的每條線段是否都在多邊形內(nèi)即可。設(shè)折線有m條線段,多邊形有n個頂點,則該算法的時間復(fù)雜度為O(m*n)。
判斷多邊形是否在多邊形內(nèi):
只要判斷多邊形的每條邊是否都在多邊形內(nèi)即可。判斷一個有m個頂點的多邊形是否在一個有n個頂點的多邊形內(nèi)復(fù)雜度為O(m*n)。
判斷矩形是否在多邊形內(nèi):
將矩形轉(zhuǎn)化為多邊形,然后再判斷是否在多邊形內(nèi)。
判斷圓是否在多邊形內(nèi):
只要計算圓心到多邊形的每條邊的最短距離,如果該距離大于等于圓半徑則該圓在多邊形內(nèi)。計算圓心到多邊形每條邊最短距離的算法在后文闡述。
判斷點是否在圓內(nèi):
計算圓心到該點的距離,如果小于等于半徑則該點在圓內(nèi)。
判斷線段、折線、矩形、多邊形是否在圓內(nèi):
因為圓是凸集,所以只要判斷是否每個頂點都在圓內(nèi)即可。
判斷圓是否在圓內(nèi):
設(shè)兩圓為O1,O2,半徑分別為r1, r2,要判斷O2是否在O1內(nèi)。先比較r1,r2的大小,如果r1<r2則O2不可能在O1內(nèi);否則如果兩圓心的距離大于r1 - r2 ,則O2不在O1內(nèi);否則O2在O1內(nèi)。
計算點到線段的最近點:
如果該線段平行于X軸(Y軸),則過點point作該線段所在直線的垂線,垂足很容易求得,然后計算出垂足,如果垂足在線段上則返回垂足,否則返回離垂足近的端點;如果該線段不平行于X軸也不平行于Y軸,則斜率存在且不為0。設(shè)線段的兩端點為pt1和pt2,斜率為:k = ( pt2.y - pt1. y ) / (pt2.x - pt1.x );該直線方程為:y = k* ( x - pt1.x) + pt1.y。其垂線的斜率為 - 1 / k,垂線方程為:y = (-1/k) * (x - point.x) + point.y 。
聯(lián)立兩直線方程解得:x = ( k^2 * pt1.x + k * (point.y - pt1.y ) + point.x ) / ( k^2 + 1) ,y = k * ( x - pt1.x) + pt1.y;然后再判斷垂足是否在線段上,如果在線段上則返回垂足;如果不在則計算兩端點到垂足的距離,選擇距離垂足較近的端點返回。
計算點到折線、矩形、多邊形的最近點:
只要分別計算點到每條線段的最近點,記錄最近距離,取其中最近距離最小的點即可。
計算點到圓的最近距離及交點坐標:
如果該點在圓心,因為圓心到圓周任一點的距離相等,返回UNDEFINED。
連接點P和圓心O,如果PO平行于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,這時直線PO斜率為k = ( P.y - O.y )/ ( P.x - O.x )。直線PO的方程為:y = k * ( x - P.x) + P.y。設(shè)圓方程為:(x - O.x ) ^2 + ( y - O.y ) ^2 = r ^2,聯(lián)立兩方程組可以解出直線PO和圓的交點,取其中離P點較近的交點即可。
計算兩條共線的線段的交點:
對于兩條共線的線段,它們之間的位置關(guān)系有下圖所示的幾種情況。圖(a)中兩條線段沒有交點;圖 (b) 和 (d) 中兩條線段有無窮焦點;圖 (c) 中兩條線段有一個交點。設(shè)line1是兩條線段中較長的一條,line2是較短的一條,如果line1包含了line2的兩個端點,則是圖(d)的情況,兩線段有無窮交點;如果line1只包含line2的一個端點,那么如果line1的某個端點等于被line1包含的line2的那個端點,則是圖(c)的情況,這時兩線段只有一個交點,否則就是圖(b)的情況,兩線段也是有無窮的交點;如果line1不包含line2的任何端點,則是圖(a)的情況,這時兩線段沒有交點。

計算線段或直線與線段的交點:
設(shè)一條線段為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共線,假如L1是直線的話他們有無窮的交點,假如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共線,假如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) 如果K1 = K2
i. 如果Q1在L0上,則說明L0和L1共線,假如L1是直線的話有無窮交點,假如L1是線段的話可用"計算兩條共線線段的交點"的算法求他們的交點(該方法在前文已討論過);
ii. 如果Q1不在L0上,則說明L0和L1平行,他們沒有交點。
c) 聯(lián)立兩直線的方程組可以解出交點來
這個算法并不復(fù)雜,但是要分情況討論清楚,尤其是當兩條線段共線的情況需要單獨考慮,所以在前文將求兩條共線線段的算法單獨寫出來。另外,一開始就先利用矢量叉乘判斷線段與線段(或直線)是否相交,如果結(jié)果是相交,那么在后面就可以將線段全部看作直線來考慮。需要注意的是,我們可以將直線或線段方程改寫為ax+by+c=0的形式,這樣一來上述過程的部分步驟可以合并,縮短了代碼長度,但是由于先要求出參數(shù),這種算法將花費更多的時間。
求線段或直線與折線、矩形、多邊形的交點:
分別求與每條邊的交點即可。
求線段或直線與圓的交點:
設(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平行于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)證明了凸包算法的時間復(fù)雜度下界是O(n*logn),但是當凸包的頂點數(shù)h也被考慮進去的話,Krikpatrick和Seidel的剪枝搜索算法可以達到O(n*logh),在漸進意義下達到最優(yōu)。最常用的凸包算法是Graham掃描法和Jarvis步進法。本文只簡單介紹一下Graham掃描法,其正確性的證明和Jarvis步進法的過程大家可以參考《算法導(dǎo)論》。
對于一個有三個或以上點的點集Q,Graham掃描法的過程如下:
令p0為Q中Y-X坐標排序下最小的點
設(shè)<p1,p2,...pm>為對其余點按以p0為中心的極角逆時針排序所得的點集(如果有多個點有相同的極角,除了距p0最遠的點外全部移除
壓p0進棧S
壓p1進棧S
壓p2進棧S
for i ← 3 to m
do while 由S的棧頂元素的下一個元素、S的棧頂元素以及pi構(gòu)成的折線段不拐向左側(cè)
對S彈棧
壓pi進棧S
return S;
此過程執(zhí)行后,棧S由底至頂?shù)脑鼐褪荙的凸包頂點按逆時針排列的點序列。需要注意的是,我們對點按極角逆時針排序時,并不需要真正求出極角,只需要求出任意兩點的次序就可以了。而這個步驟可以用前述的矢量叉積性質(zhì)實現(xiàn)。
四、結(jié)語
盡管人類對幾何學的研究從古代起便沒有中斷過,但是具體到借助計算機來解決幾何問題的研究,還只是停留在一個初級階段,無論從應(yīng)用領(lǐng)域還是發(fā)展前景來看,計算幾何學都值得我們認真學習、加以運用,希望這篇文章能帶你走進這個豐富多彩的世界。
posted @
2009-10-06 14:38 wyiu 閱讀(179) |
評論 (0) |
編輯 收藏
500MS弄到79MS,求多邊形面積,有可能是凹,用叉積就可以了,不能存所有點,會MLE
一次讀入比單個讀入快
#include<iostream>
#include<math.h>
using namespace std;

struct point
{__int64 x,y;};

__int64 dx[]=
{0,-1,0,1,-1,0,1,-1,0,1};

__int64 dy[]=
{0,-1,-1,-1,0,0,0,1,1,1};

//__int64 xmult(point p0,point p1,point p2){
// return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
//}
char path[1000010];
int main()


{
int t;
point p0,p1,p2;
__int64 ans;
scanf("%d",&t);
while(t--)

{
p0.x=0;
p0.y=0;
p2=p1=p0;
ans=0;
scanf(" %s",path);
int i=0;
for(i=0;path[i]!='5';i++)

{
p1=p2;
p2.x+=dx[path[i]-'0'];
p2.y+=dy[path[i]-'0'];
ans+=(p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
if(ans<0) ans*=-1;
printf("%I64d",(ans>>1));
if(((ans>>1)<<1)==ans) printf("\n");
else printf(".5\n");
}
return 0;
}

posted @
2009-10-06 13:12 wyiu 閱讀(138) |
評論 (0) |
編輯 收藏
//枚舉4條邊界的點a,然后求線段out(trea,a)與其他wall相交的最小數(shù)
//代碼中的max改為min有助理解
#include<iostream>
using namespace std;
#define MAXN 30
#define eps 1e-8
#define zero(x) (((x)>0?(x):-(x))<eps)


struct point
{double x,y;};

struct line
{point a,b;};

line wall[MAXN+1];
int n;
point trea;
//計算cross product (P1-P0)x(P2-P0)

double xmult(point p1,point p2,point p0)
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}

//判兩點在線段異側(cè),點在線段上返回0

int opposite_side(point p1,point p2,line l)
{
return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)<-eps;
}
//判兩線段相交,不包括端點和部分重合

int intersect_ex(line u,line v)
{
return opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);
}
int main()


{
int i;
scanf("%d",&n);
for(i=0;i<n;i++)
scanf("%lf%lf%lf%lf",&(wall[i].a.x),&(wall[i].a.y),&(wall[i].b.x),&(wall[i].b.y));
scanf("%lf%lf",&(trea.x),&(trea.y));
line out;
out.a=trea;
int max=n+1;
double x,y;
for(x=0,out.b.y=100;x<=100;x+=0.01)

{
out.b.x=x;
int t=0;
for(i=0;i<n;i++)
if(intersect_ex(out,wall[i]))
t++;
if(t<max) max=t;
}
for(x=0,out.b.y=0;x<=100;x+=0.01)

{
out.b.x=x;
int t=0;
for(i=0;i<n;i++)
if(intersect_ex(out,wall[i]))
t++;
if(t<max) max=t;
}
for(y=0,out.b.x=0;y<=100;y+=0.01)

{
out.b.y=y;
int t=0;
for(i=0;i<n;i++)
if(intersect_ex(out,wall[i]))
t++;
if(t<max) max=t;
}
for(y=0,out.b.x=100;y<=100;y+=0.01)

{
out.b.y=y;
int t=0;
for(i=0;i<n;i++)
if(intersect_ex(out,wall[i]))
t++;
if(t<max) max=t;
}
printf("Number of doors = %d\n",max+1);
return 0;
}
posted @
2009-10-05 15:05 wyiu 閱讀(134) |
評論 (0) |
編輯 收藏
#include<stdio.h>
#include<math.h>
#include<memory.h>
//using namespace std;
#define eps 1e-8
#define MAXN 50
#define PI 3.1415926535898

struct point
{ double x,y;};

point planet[MAXN+1];
bool visit[MAXN+1];
int ans[MAXN+1];


double xmult(point p0,point p1,point p2)
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}

double dmult(point p0,point p1,point p2)
{
return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y);
}


double distance(point p1,point p2)
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}

int main()


{
int t,indice,i,j,k,n;
double max;
scanf("%d",&t);
while(t--)

{
scanf("%d",&n);
for(i=1;i<=n;i++)
scanf("%d%lf%lf",&indice,&(planet[i].x),&(planet[i].y));

memset(visit,0,sizeof(visit));

int next=1;

for(i=2;i<=n;i++)
if(planet[next].y>planet[i].y)
next=i;
else if(planet[next].y == planet[i].y)
if(planet[next].x>planet[i].x)
next=i;
point p1,p0;
p1.x=0;
p1.y=planet[next].y;
p0.x=planet[next].x;
p0.y=planet[next].y;
k=0;
double thta,max;
for(j=1;j<=n;j++)

{
visit[next]=true;
ans[k++]=next;
max=2*PI;
for(i=1;i<=n ; i++)

{
if(visit[i]) continue;
double xmul=xmult(p0,p1,planet[i]);
double dmul=dmult(p0,p1,planet[i]);
if( xmul<=0)

{
if(dmul==0 )
thta=3*PI/2.0;
else
if(dmul<0 ) thta=PI+atan(xmul/dmul);
else thta=2*PI+atan(xmul/dmul);

}
if( thta<max ||
(fabs(thta-max)<eps && fabs(distance(p0,planet[i])-distance(p0,planet[next]))<eps)
)

{ next=i;max=thta;}
}
p1.x=p0.x;
p1.y=p0.y;
p0.x=planet[next].x;
p0.y=planet[next].y;
}
printf("%d",k);
for(i=0;i<k;i++)
printf(" %d",ans[i]);
printf("\n");
}
return 0;
}
posted @
2009-10-05 12:12 wyiu 閱讀(111) |
評論 (0) |
編輯 收藏
// 同樣是暴力,不同的暴力,一個TLE,一個5**MS
//判
線段相交
#include<iostream>
#include <math.h>
using namespace std;
#define MAXN 100000
#define eps 1e-8
#define zero(x) (((x)>0?(x):-(x))<eps)


struct point
{double x,y;};

struct line
{ point a,b;};
line sti[MAXN+1];
bool under[MAXN+1];

//計算cross product (P1-P0)x(P2-P0)

double xmult(point p1,point p2,point p0)
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}


int dots_inline(point p1,point p2,point p3)
{
return zero(xmult(p1,p2,p3));
}

int same_side(point p1,point p2,line l)
{
return xmult(l.a,p1,l.b)*xmult(l.a,p2,l.b)>eps;
}

int dot_online_in(point p,line l)
{
return zero(xmult(p,l.a,l.b))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&(l.a.y-p.y)*(l.b.y-p.y)<eps;
}



int intersect_in(line u,line v)
{
if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))
return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);
return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);
}


int main()


{
int n,i,j;
while(scanf("%d",&n)!=EOF && n)

{
memset(under,0,sizeof(under));
for(i=1;i<=n;i++)
scanf("%lf%lf%lf%lf",&(sti[i].a.x),&(sti[i].a.y),&(sti[i].b.x),&(sti[i].b.y));
for(i=1;i<n;i++)
for(j=i+1;j<=n;j++)

{
if(intersect_in(sti[i],sti[j]))

{
under[i]=true;
break;
}
}
printf("Top sticks:");
for(i=1;i<n;i++)
if(!under[i])
printf(" %d,",i);
printf(" %d.\n",n);
}
return 0;
}
posted @
2009-10-04 21:39 wyiu 閱讀(214) |
評論 (0) |
編輯 收藏
一.叉積
設(shè) a(x1,y1), b(x2,y2)
二維:a x b=x1*y2-x2*y1
設(shè)p0(x0,y0), p1(x1,y1) ,p2(x2,y2), p3(x3,y3)
< p0p1>= (p1-p0) =(x1-x0,y1-y0);
<p2p3>= (p3-p2) =(x3-x2,y3-y2);
<p0p1> x <p2p3> =(p1-p0) x (p3-p2) = (x1-x0)*(y3-y2)- (x3-x2)* (y1-y0);
(p1-p0) x (p3-p2)結(jié)果的意義:
正: <p0p1>在<p2,p3>順時針(0,pi)內(nèi)
負: <p0p1>在<p2,p3>逆時針(0,pi)內(nèi)
0 : <p0p1>, <p2p3>共線,夾角為0或pi

double xmult(point p0,point p1,point p2)
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}

可用于判p2與<p0 p1>的關(guān)系

double xmult(point p0,point p1,point p2,point p3)
{
return (p1.x-p0.x)*(p3.y-p2.y)-(p3.x-p2.x)*(p1.y-p0.y);
}
可用于判<p0p1>, <p2p3>位置關(guān)系

double xmult(line u,line v)
{
return (u.b.x-u.a.x)*(v.b.y-v.a.y)-(v.b.x-v.a.x)*(u.b.y-u.a.y);
}
//計算兩直線交點,注意事先判斷直線是否平行!

point intersection(line u,line v)
{
point ret=u.a;
double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))
/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));
ret.x+=(u.b.x-u.a.x)*t;
ret.y+=(u.b.y-u.a.y)*t;
return ret;
}
int relat_and_intersection(line u,line v,point &ret)


{

/**//*
LINE 重合
NONE 平行
INTERSECT 相交并返回交點
*/
double a1,b1,c1,a2,b2,c2;
//系數(shù)
a1=u.a.y-u.b.y;
b1=-(u.a.x-u.b.x);
c1=u.a.x*u.b.y-u.a.y*u.b.x;

a2=v.a.y-v.b.y;
b2=-(v.a.x-v.b.x);
c2=v.a.x*v.b.y-v.a.y*v.b.x;

if(fabs(a1*b2-a2*b1)<eps)

{
if( fabs(a1*c2-a2*c1)<eps && fabs(b1*c2-b2*c1)<eps )
return LINE;
else return NONE;
}
else

{
point ret;
ret.x=(b1*c2-b2*c1)/(a1*b2-a2*b1);
ret.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);
return INTERSECT;
}
}
posted @
2009-10-04 18:38 wyiu 閱讀(254) |
評論 (0) |
編輯 收藏
//學習判定兩直線關(guān)系:平行 重合 相交
//求交點
#include<iostream>
#include <math.h>
using namespace std;
#define eps 1e-8
#define zero(x) (((x)>0?(x):-(x))<eps)

struct point
{double x,y;};

struct line
{point a,b;};


double xmult(line u,line v)
{
return (u.b.x-u.a.x)*(v.b.y-v.a.y)-(v.b.x-v.a.x)*(u.b.y-u.a.y);
}
//計算兩直線交點,注意事先判斷直線是否平行!

point intersection(line u,line v)
{
point ret=u.a;
double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))
/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));
ret.x+=(u.b.x-u.a.x)*t;
ret.y+=(u.b.y-u.a.y)*t;
return ret;
}

int main()


{
line u,v;
int n;
scanf("%d",&n);
printf("INTERSECTING LINES OUTPUT\n");
while(n--)

{
scanf("%lf%lf%lf%lf",&(u.a.x),&(u.a.y),&(u.b.x),&(u.b.y));
scanf("%lf%lf%lf%lf",&(v.a.x),&(v.a.y),&(v.b.x),&(v.b.y));
double a1,b1,c1,a2,b2,c2;
//xishu
a1=u.a.y-u.b.y;
b1=-(u.a.x-u.b.x);
c1=u.a.x*u.b.y-u.a.y*u.b.x;

a2=v.a.y-v.b.y;
b2=-(v.a.x-v.b.x);
c2=v.a.x*v.b.y-v.a.y*v.b.x;
//if( fabs(xmult(u,v))<eps)
if(fabs(a1*b2-a2*b1)<eps)

{
if( fabs(a1*c2-a2*c1)<eps && fabs(b1*c2-b2*c1)<eps )
printf("LINE\n");
else printf("NONE\n");
}
else

{
point res;
//res=intersection(u,v);
res.x=(b1*c2-b2*c1)/(a1*b2-a2*b1);
res.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);
printf("POINT %.2lf %.2lf\n",res.x,res.y);
}
}
printf("END OF OUTPUT\n");
return 0;
}
posted @
2009-10-04 18:15 wyiu 閱讀(149) |
評論 (0) |
編輯 收藏