筆者極端厭惡計算幾何,以致於直到今天以前都只打過幾個計算幾何的模板~~~~~算法
不過鑑於今年是18年,因此感受SD頗有可能考計算幾何(18年是什麼理由啊喂)數組
因而滾過來整理計算幾何的資料了......ide
《計算幾何——算法與應用》學習
Mark·Overmars 著 鄧俊輝 譯優化
《計算幾何導論》搜索引擎
[美] F·P·普霍帕拉塔 M·I·沙莫斯 著spa
——好像是一些沒什麼用的東西呢3d
《算法導論——第33章》code
[美]Tomas H.Cormen Charles E.Leiserson Ronald L.Rivest Clifford Stein 著orm
——據說是能夠提高文章逼格的東西呢
本文會盡可能(大概吧)把上述文獻中對OI有用的部分扒拉出來
《藍書》
劉汝佳 著
《紫書》
劉汝佳 著
某些大佬的blog——%%%%%%%%%%%%%%%%%%%%%%%%%%%
BZOJ——儘管他很爆炸
POJ——儘管全是英文
HDU——我好像沒怎麼用過
luogu——很良心的OJ吧,大概......
LOJ——幾乎沒用過,應該也用不到吧
UOJ——不用搜索引擎就能翻出須要的題——若是須要的題真的湊巧在這個OJ裏的話
(廢話真多)
計算幾何最重要的是(拼命調)關注點和向量
計算幾何的研究對象每每是基於歐幾里得空間的點集的,
由於,當咱們存在一個座標系時
一個點能夠有一個座標,因而與原點構成了向量——實際上,在許多大佬的計算幾何模板中都將點和向量宏定義爲一個東西,如劉汝佳的《藍書》
兩個點能夠構成一個線段或者直線,直線能夠定義半平面交
多個點能夠構成多邊形,或者凸包——事實上凸包和半平面交存在着諸多類似性
多邊形能夠各類剖分
並且,相比於用其餘元素描述上述內容而言,座標系中的點是更加適用於計算機語言的,由於一個點表示成有序數對以後,與任何數對(pair)無異
固然以上胡扯也能夠扯到擴展到更高的維度上去
因而約定有以下不經常使用的記法:
$E^d$表示d維歐幾里得空間,(在下文中,d通常大於等於2且不超過3,——因爲筆者水平而有的限制)
在d維歐幾里得空間中,有以下基本對象:
一個d元組$(x_1,x_2,...,x_d)$,他能夠表示一個點P,也能夠表示一個起點爲$E^d$的原點O,終點爲點P的一個有d個份量的向量
給定$E^d$中兩個不一樣的點$P_1$,$P_1$的線性組合
$$aP_1+(1-a)P_2(a∈R)$$
能夠表示$E^d$中的一條直線(若是把P看作向量的話,上式是高中數學中常見的向量基本定理推論)
給定$E^d$中三個不一樣的點$P_1$,$P_2$,$P_3$的線性組合
$$aP_1+bP_2+(1-a-b)P_3$$
能夠表示$E^d$中的一個平面
相似的上面關於線面的定義也能夠推廣的更高的維度而成爲一種叫作線性簇的東西
(沒有用)
給上文中的直線表達式中加一個對a的範圍限制,能夠獲得一條線段
如,加入限制(0≤a≤1)能夠獲得一個以$P_1$,$P_2$位端點的線段,能夠表示爲$\overline{P_1P_2}$或$\overline{P_2P_1}$
若$E^d$中的點集D中的任何兩點$P_1,P_2$,$\overline{P_1P_2}$中的全部點屬於D,則稱D是$E^d$中的一個凸集
凸集的交是凸集
直線中的任意線段是凸集
平面內的任意凸多邊形是凸集
空間中的任意凸多面體是凸集
關於更高維的凸集的相關性質,可移步《線性規劃》——反正筆者看了一點就直接棄療了的說
$E^d$中點集S的凸包是$E^d$中包含S的最小凸集,
凸殼是凸包的邊界
在$E^2$中多邊形定義爲線段的一個有限集合使得每一個端點恰爲兩條邊所共有,且沒有邊的子集具備該性質(多邊形僅指這一個邊界部分)
注意:
是多邊形
而
不是
多邊形的其餘內容沒有什麼值得強調的
若不相鄰的邊無公共點,則多邊形是簡單的
如上上圖中的多邊形就不是簡單的
簡單多邊形把平面劃成兩個不相交的區域
稱爲爲多邊形的內部和外部
若簡單多邊形P的內部是個凸集,則此簡單多邊形是凸多邊形
若存在點z在簡單多邊形P的內部或邊上,知足對於P的全部點p,線段$\overline{zp}$屬於P的內部或屬於P邊,則說明多邊形P是星形的
(所以每一個凸多邊形都是星形的)
然而:
如上圖就不是星形多邊形
在星形多邊形P中,全部z構成的集合爲P的核(he?hu?)容易證實,核是一個凸集
(能夠用凸集的定義來證)
這一部分將在$E^2$內
以向量的計算爲基礎
討論點與線段的關係,線段與線段的關係
將介紹如何定義一個點,如何定義一個向量
如何實現向量的基本計算
如何用點和向量表示直線、線段,進而用點和向量實現對直線線段的操做
這大概是計算幾何中最簡單的一部分,將會粘貼一些模板
向量的基本運算做爲處理計算幾何點線的基礎出現
const double eps=1e-10; int dcmp(double x){ if(fabs(x)<eps)return 0; else return x<0?-1:1; }//精度
struct Point{ double x,y; Point(double x=0,double y=0):x(x),y(y){ }; };//定義點 typedef Point Vector;//定義向量
Vector operator + (Point A,Point B){ return Vector(A.x+B.x,A.y+B.y); }//定義向量加符合三角形定則
Vector operator - (Point A,Point B){ return Vector(A.x-B.x,A.y-B.y); }//定義向量減符合三角形定則
Vector operator * (Vector A,double p){ return Vector(A.x*p,A.y*p); }//定義向量數乘 Vector operator * (double p,Vector A){ return Vector(A.x*p,A.y*p); }//定義向量數乘
(直接分別定義左右乘了,並不會其餘更高端的方式)
Vector operator / (Vector A,double p){ return Vector(A.x/p,A.y/p); }//定義向量數除
bool operator == (const Point& a,const Point& b){ return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; }//定義向量相等爲徹底相同
即從x軸正半軸旋轉到該向量方向所需的弧度,(單位:弧度)
double Polar_Angle (Vector A){ return atan2(A.y,A.x); }//計算向量極角
將向量p(x,y)繞起點逆時針旋轉弧度a,獲得$p'(xcosa-ysina,xsina+ycosa)$
若是你已經學習了《數學必修》的內容,則你應該對點積十分了解
double Dot(Vector A,Vector B){ return A.x*B.x+A.y*B.y; }//計算向量點積
座標對應相乘再相加,點積的結果是一個實數,可認爲是模長乘夾角cos
若是你有必定的物理或數學水平,則你應該明白二維向量的叉積結果應該是三維空間中的一個垂直於這個兩個向量向量,方向按照右手法則斷定
然而由於使用範圍的限制,這裏對叉積的定義十分淺顯:
考慮向量$p_1(x_1,y_1),p_2(x_1,y_1)$的叉積定義爲把兩個向量起點移動到同一處後獲得的三點圍成的三角形的有向面積的兩倍
或者定義爲下述矩陣的行列式的值:
若其爲正,則相對於原點來講$p_1$位於$p_2$的順時針方向,若其爲0,則兩向量共線
double Cross(Vector A,Vector B){ return A.x*B.y-A.y*B.x; }//計算向量叉積
交叉相乘再相減,可看作模長乘夾角sin
注意:叉積沒有交換律
因而有了一個簡單的基於叉積的計算三點圍成三角形面積的算法:
double Area(Point A,Point B,Point C){ return fabs(Cross(B-A,C-A)/2); }//計算abc三點三角形的數值面積
因而有了一個簡單的基於點積的計算向量長的算法:
double Length(Vector A){ return sqrt(Dot(A,A)); }//計算向量模長
因而有了一個簡單的基於點積的計算向量夾角的算法:
double Angle(Vector A,Vector B){ return acos(Dot(A,B)/Lenth(A)/Lenth(B)); }//計算向量夾角(有序的)
線段取其兩端點,定義爲無序的一對點
而直線則取其中一點和直線的方向向量,定義爲一個點和一個向量的二元組
(也能夠看作兩個向量,也能夠看作兩個點)
這一部分的存在應用價值(嗎?)
討論在點$p_1$處,連續線段$\overline{p_0p_1}$,$\overline{p_1p_2}$的旋轉方向,即斷定給定角$∠p_0p_1p_2$的轉向,而咱們能夠經過計算$\overrightarrow{p_0p_1}$$\overrightarrow{p_0p_2}$叉積結果,來避免直接對角的計算,若其爲正,則說明後一條線段向逆時針偏折;若其爲負,則說明後一條線段向順時針偏折;若其爲零,則兩條線段共線
int Direction(Vector p0,Vector p1,Vector p2){ double num=Cross(Vector(p1-p0),Vector(p2-p0)); return dcmp(num); }//判斷p0->p1,p1->p2的偏轉方向,1爲逆時針,-1爲順時針,0爲共線
在以前,咱們定義的線性組合$a\vec{A}+(1-a)\vec{B}(a∈R)$爲直線,
簡單地,咱們定義直線的方向向量爲$\vec{v}=\vec{A}-\vec{B}$,
化簡,因而有了$B+a\vec{v}$
方向向量界定了直線的方向,向量B的終點B界定了直線在笛卡爾座標系中的位置,使用規範的符號,因而有了參數方程:
$$P+t\vec v$$
其中t的不一樣取值表明了直線上的不一樣點
因而,設$l_1:P+t\vec v $,$l_2:Q+t\vec w $,聯立兩參數方程可解得交點在$l_1$的參數$t_1$,在$l_2$的參數$t_2$,若設$\vec u =P-Q$則有關於$t_1,t_2$的公式以下:
$$t_1={cross(\vec w,\vec u)\over cross(\vec v,\vec w)},t_2={cross(\vec v,\vec u)\over cross(\vec v,\vec w)}$$
Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){ Vector u=P-Q; double t=Cross(w,u)/Cross(v,w); return P+v*t; }//用參數方程求取直線交點——若是他存在的話
實際上,運用簡單的高中解析幾何知識就足以解決此問題了吧
這一部分有應用的價值(嗎)
直接叉積求解便可
double DistanceToLine(Point P,Point A,Point B){ return abs(Cross(B-A,P-A)/Length(B-A)); }//P到直線AB的距離
須要使人十分不愉快的分類討論:
設P在直線AB上的投影點爲Q,
\
1:Q在線段外,據A較近,須要求解PA的長
2:Q在線段AB上,直接使用點到直線距離
3:Q在線段外,據B較近,須要求解PB的長
斷定能夠經過點積來實現
double DistanceToSigment(Point P,Point A,Point B){ if(A==B)return Length(P-A); if(dcmp(Dot(P-A,B-A))<0) return Length(P-A); if(dcmp(Dot(P-B,A-B))<0) return Length(P-B); return DistanceToLine(P,A,B); }//P到線段AB的距離
習題:
POJ P2318
POJ P1269
UVa11796(UVA的小題目挺有意思的,適合心情好的時候)
因而咱們很愉快地從對點線的討論正式轉入圖形了
這一部分將討論對多邊形問題的解決,
進而研究凸包凸殼半平面交等基於點線對象
而後將給出圓的定義
多邊形是在二維幾何中常常研究的圖形,
多邊形取其全部頂點和他們的逆時針鏈接順序而定義爲一個有序點集
具體的定義已經在本文的第一部分說起,因此不做贅述,
接下來將介紹簡單多邊形圍成有向面積的求法
點與多邊形關係的斷定
考慮凸多邊形的有向面積?
能夠從一個點出發把凸多邊形劃成n-2個三角形,求面積再求和
這個方法對非凸的簡單多邊形一樣適用——由於三角形面積的有向性使得在多邊形外部的部分被抵消掉
具體的作法是,把多邊形的n個頂點按照鏈接順序從0到n-1編號,從$P_0$依次向$P_1$到$P_{n-1}$連向量$\vec{v}$,計算$\sum_{i=1}^{n-2}{cross(\vec{v_i},\vec{v_{i+1}})\over 2}$
(深色部分正負相抵)
double PolygonArea(Point *p,int n){ double area=0; int i; for(i=1;i<n-1;i++) area+=Cross(p[i]-p[0],p[i+1]-p[0]); return area/2; }//計算多邊形的有向面積,也許在結果上加個||更有用些??
一下算法適用於全部點按鏈接順序逆時針排列的簡單多邊形
如何斷定p點是否在多邊形Poly內部呢
粗略地想,從點p向右引一條平行於x軸的射線
因爲射線無限長而多邊形內部範圍有限,
因此射線無限遠處的盡頭所在的區域必定是多邊形外部
因爲多邊形每一條邊兩側必定分屬多邊形的內外
因此從射線遠離p的一側向p走,
第一次通過一條邊將到達多邊形內
第二次通過一條邊將到達多邊形外
以此類推通過——
若該射線奇數次穿過邊,則點p在多邊形內,反之亦然
一些細節,
若是射線的某段與某條邊重合,則這條邊不該計入次數
若是射線穿過某端點,應該如何斷定呢
比較好的方法是定義邊的方向爲從編號小的點到編號大的點
這樣之後再定義
從下往上穿過射線的邊包含起點不包含終點
從上往下穿過射線的邊包含終點不包含起點
這樣若穿過的端點所在的兩邊同向則只被計算一次
若穿過的端點所在的邊反向則要麼一次都不計算,要麼直接算兩次
int isPointInPolygon(Point p,Point poly[],int n){ int wn=0,i,k,d1,d2; for(i=0;i<n;i++){ if(!dcmp(DistanceToSigment(p,poly[i],poly[(i+1)%n]))) return -1;//在邊界上 k=dcmp(Cross(poly[(i+1)%n]-poly[i],p-poly[i])); d1=dcmp(poly[i].y-p.y); d2=dcmp(poly[(i+1)%n].y-p.y); if(k>0&&d1<=0&&d2>0)wn++;//點在左,下上穿 if(k<0&&d2<=0&&d1>0)wn++;//點在右,上下穿 } return wn&1;//T:內部;F:外部 }//斷定點是否在多邊形內部
凸多邊形的對踵點是凸多邊形上存在的不僅一對的點對,
每對對踵點$(p_1,p_1)$知足都是凸多邊形的頂點且存在過$p_1$的某條凸多邊形的切線與某條過$p_2$的凸多邊形切線平行
對踵點個數不超過(3n/2)
旋轉卡殼算法是一種經過處理凸多邊形的對踵點來解決一系列凸多邊形問題的算法
筆者將旋轉卡殼算法從多邊形中剝離出來(主要是由於他內容比較多)
這裏給出用旋轉卡殼枚舉全部對踵點的方法:
先找到y最小的點$p=p_{ymin}$,和y最大的點$q=p_{ymax}$,他們是一對對踵點,
從這對對踵點分別向逆時針方向引出平行於x軸的射線,
按相同的角速度逆時針轉動這對射線,直到其中一條穿過多邊形的下一端點$p_{next}$,{p_{next}}將代替他所在的射線的原端點成爲射線的新端點,並和另外一射線的端點構成新的對踵點,
繼續旋轉新的這對射線以得到新的對踵點
當咱們枚舉了全部角度的平行線以後,天然也就枚舉出了全部對踵點
枚舉角度的效率是INF的(十分高效)
然而並非全部角度的平行切線都切與不一樣的對踵點,因此在轉動的過程當中有許多角度能夠直接略過
具體地說,同一角速度轉動這一說法太過玄學,具體實現應該是,
把上述過程按照兩個射線的頂點是否有其一被替代來劃分爲許多階段,
能夠經過叉積來斷定$\overrightarrow {pp_{next}}$和$\overrightarrow{q_{next}q}$的極角大小,
進而肯定p與q被替代的前後順序
而後直接在兩個階段之間跳躍,由於這期間沒有新的對踵點產生
凸多邊形的直徑是凸多邊形邊上全部點中距離最遠的一對點的距離,
他顯然是某對對踵點的距離,枚舉全部對踵點便可
效率$O(n)$(對每一對對踵點O(1))
double RotatingCaliper_diameter(Point poly[],int n){ int p=0,q=n-1,fl,i; double ret=0; for(i=0;i<n;i++){ if(poly[i].y>poly[q].y) q=i; if(poly[i].y<poly[p].y) p=i; } Point ymin=poly[p],ymax=poly[q]; for(i=0;i<n*3;i++){ if(ret<Length(poly[p%n]-poly[q%n])) ret=Length(poly[p%n]-poly[q%n]); fl=dcmp(Cross(poly[(p+1)%n]-poly[p%n],poly[q%n]-poly[(q+1)%n])); if(!fl){ if(ret<Length(poly[p%n]-poly[(q+1)%n])) ret=Length(poly[p%n]-poly[(q+1)%n]); if(ret<Length(poly[q%n]-poly[(p+1)%n])) ret=Length(poly[q%n]-poly[(p+1)%n]); p++,q++; } else{ if(fl>0) p++; else q++; } } return ret; }//用旋轉卡殼求解凸多邊形直徑
習題:
POJ P2187——須要求凸包,而後求凸包的直徑的平方輸出
凸多邊形的寬度是凸多邊形的全部平行切線對的距離的最小值
對任意一對平行切線$l_1^a,l_2^a$,必有一對平行切線$l_1^b,l_2^b$,其中至少一條與某一整條邊重合
知足$dis(l_1^b,l_2^b)<=dis(l_1^a,l_2^a)$
因而枚舉此類平行切線便可
因爲,旋轉卡殼的過程正是由出現此類平行切線的時間點來劃分的,因此能夠枚舉全部此類平行切線
效率$O(n)$(對每一對對踵點O(1))
double RotatingCaliper_width(Point poly[],int n){ int p=0,q=n-1,fl,i; double ret; for(i=0;i<n;i++){ if(poly[i].y>poly[q].y) q=i; if(poly[i].y<poly[p].y) p=i; } ret=Length(poly[p]-poly[q]); for(i=0;i<n*3;i++){ fl=dcmp(Cross(poly[(p+1)%n]-poly[p%n],poly[q%n]-poly[(q+1)%n])); if(!fl){ if(ret>DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n])) ret=DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n]); p++,q++; } else{ if(fl>0){ if(ret>DistanceToLine(poly[q%n],poly[p%n],poly[(p+1)%n])) ret=DistanceToLine(poly[q%n],poly[p%n],poly[(p+1)%n]); p++; } else{ if(ret>DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n])) ret=DistanceToLine(poly[p%n],poly[q%n],poly[(q+1)%n]); q++; } } } return ret; }//用旋轉卡殼求解凸多邊形寬度
(沒找題,正好yhzq大神作的某次比賽有這個題,正好用那個題實驗了板子的正確性2333)
這樣的矩形必定與凸多邊形至少一邊發生重合,因而這一條重合的邊以及這條邊的對邊能夠經過旋轉卡殼來枚舉,
在枚舉出這麼一組有重合平行邊以後,如何枚舉另一對平行邊呢,
在咱們枚舉出第一對有重合平行邊並找到其對應的另外一對邊以後(這個另外一對邊方向與有重合平行邊垂直,因此這對邊其實能夠存成多邊形上的一對點)
若是咱們旋轉這組有重合平行邊獲得另外一組有重合平行邊的話,
新的一組的對應另外一對邊能夠由舊的一組經過向相同的方向旋轉得來,
斷定是否完成旋轉的方法能夠是向量叉積
效率$O(n)$(對每一對對踵點O(1),另外一對平行線的旋轉也是單向的因而也是O(n))
其實求最小周長外接矩形也是同理
void RC(Poi poly[],int n){ int p=0,q=n-1,l=0,r=n-1; int fl,i,j; Vec nm,dr; for(i=0;i<n;i++){ if(poly[i].y<poly[p].y) p=i; if(poly[i].y>poly[q].y) q=i; if(poly[i].x<poly[l].x) l=i; if(poly[i].x>poly[r].x) r=i; } an[0]=intersect_p(poly[p],Vec(1,0),poly[l],Vec(0,1)),an[1]=intersect_p(poly[p],Vec(1,0),poly[r],Vec(0,1)); an[2]=intersect_p(poly[r],Vec(0,1),poly[q],Vec(1,0)),an[3]=intersect_p(poly[l],Vec(0,1),poly[q],Vec(1,0)); ans=fabs(Area(an[0],an[1],an[2])); for(i=0;i<n*3;i++){ fl=dcmp(Cross(poly[(p+1)%n]-poly[p%n],poly[q%n]-poly[(q+1)%n])); if(!fl){ dr=poly[(p+1)%n]-poly[p%n],dr=dr/Len(dr); nm=Vec(dr.y,-dr.x); while(dcmp(Cross(poly[(l+1)%n]-poly[l%n],nm))>0) l++; nm=Vec(0,0)-nm; while(dcmp(Cross(poly[(r+1)%n]-poly[r%n],nm))>0) r++; cha[0]=intersect_p(poly[p%n],dr,poly[l%n],nm),cha[1]=intersect_p(poly[p%n],dr,poly[r%n],nm); cha[2]=intersect_p(poly[r%n],nm,poly[q%n],dr),cha[3]=intersect_p(poly[l%n],nm,poly[q%n],dr); if(fabs(Area(cha[0],cha[1],cha[2]))<ans){ for(j=0;j<4;j++) an[j]=cha[j]; ans=fabs(Area(an[0],an[1],an[2])); } } else{ if(fl>0){ dr=poly[(p+1)%n]-poly[p%n],dr=dr/Len(dr); nm=Vec(dr.y,-dr.x); while(dcmp(Cross(poly[(l+1)%n]-poly[l%n],nm))>0) l++; nm=Vec(0,0)-nm; while(dcmp(Cross(poly[(r+1)%n]-poly[r%n],nm))>0) r++; cha[0]=intersect_p(poly[p%n],dr,poly[l%n],nm),cha[1]=intersect_p(poly[p%n],dr,poly[r%n],nm); cha[2]=intersect_p(poly[r%n],nm,poly[q%n],dr),cha[3]=intersect_p(poly[l%n],nm,poly[q%n],dr); if(fabs(Area(cha[0],cha[1],cha[2]))<ans){ for(j=0;j<4;j++) an[j]=cha[j]; ans=fabs(Area(an[0],an[1],an[2])); } p++; } else{ dr=poly[(q+1)%n]-poly[q%n],dr=dr/Len(dr); nm=Vec(-dr.y,dr.x); while(dcmp(Cross(poly[(l+1)%n]-poly[l%n],nm))>0) l++; nm=Vec(0,0)-nm; while(dcmp(Cross(poly[(r+1)%n]-poly[r%n],nm))>0) r++; cha[0]=intersect_p(poly[p%n],dr,poly[l%n],nm),cha[1]=intersect_p(poly[p%n],dr,poly[r%n],nm); cha[2]=intersect_p(poly[r%n],nm,poly[q%n],dr),cha[3]=intersect_p(poly[l%n],nm,poly[q%n],dr); if(fabs(Area(cha[0],cha[1],cha[2]))<ans){ for(j=0;j<4;j++) an[j]=cha[j]; ans=fabs(Area(an[0],an[1],an[2])); } q++; } } } }
習題:
bzoj P1185==luogu P3187 [HNOI2007]最小矩形覆蓋
(luogu沒有SPJ,卡精度;BZOJ輸出九個nan不知還能不能過)
凸包的定義前面已經給出了,並無什麼值得強調的
值得一提的是某個問題——
定義一個物品集有n個物品,每一個物品徹底由一個大小爲m的元素集中的一些元素按照某個比例融合而成
設在i物品中j元素佔比爲$x_{i,j}$
則有$\sum_{j=1}^mx_{i,j}=1(i=1...n)$
每次詢問給出一種新的物品中各元素的佔比,問他可否被原來物品集中的物品按必定比例融合來獲得
這個問題在m爲2時能夠看作一個點是否在凸包內的斷定問題
m爲3時能夠看作一個點是否在一個三維凸包內的斷定問題
。。。
好吧,咱們仍是先看看二維凸包的求法吧
求凸包,實際是求在凸包邊上的端點,或者說,是求凸殼的過程
其流程是:
先把全部點按x爲第一關鍵字,y爲第二關鍵字排序而後去重
而後將第一第二個點加入凸殼,
從第三個點開始,
當新加入多是凸殼第i個點的$p_i$時,斷定$p_{i-1}$是否須要被剔出凸殼
若連續線段$\overline{p_{i-2}p_{i-1}}$和$\overline{p_{i-1}p_i}$是向順時針偏折的,則把$p_{i-1}$剔出凸包,
而後繼續斷定$p_{i-2}$是否須要被剔出凸包......直到當咱們斷定$p_j$時,發現他不須要被剔出凸包,
而後把$p_i$放入凸殼,發現他如今多是凸殼第j+1個點(固然咯,也可能什麼都不是)
而後繼續加入凸殼的第j+2個點,重複上述過程
當咱們結束對n的處理後,咱們獲得了一個下凸殼——他的每條線都成下凸的趨勢
這樣再從n開始倒着來一遍就求出了上凸殼,合起來就是完整的凸殼啦
時間複雜度爲排序的O(nlogn).
int ConvexHull(Point*inp,int n,Point*outp){ sort(inp,inp+n); int m=-1,k,i,j; for(i=0;i<n;i++){ while(m>0&&Direction(outp[m-1],outp[m],inp[i])<=0)m--; outp[++m]=inp[i]; } k=m; for(i=n-2;i>=0;i--){ while(m>k&&Direction(outp[m-1],outp[m],inp[i])<=0)m--; outp[++m]=inp[i]; } if(m&&outp[m]==outp[0])m--; return m+1; }//假設輸入的inp數組已經去太重了,我還要寫個離散嗎(笑)
習題:
HDU1348——讀入n個點,和一個R,求n個點凸殼的周長+$2\pi R$做爲答案
通常須要排序的東西,搞成動態的話,好像都是套個堆或者平衡樹的吧
用set|平衡樹維護凸包能夠支持動態加入新點
(若是隻有刪除操做的話,能夠離線倒序變成加入)
不能同時支持刪除......
具體的說,在set分別維護上下凸殼中的點按先x再y的順序排完序後的結果,
當加入一個點後,先判斷他可能加入上凸殼仍是下凸殼,
而後若是他須要被加入凸殼的話,找到他從按先x再y的順序排完序後應該在位置,
判斷他在set中的前驅和後繼是否須要被刪除,相似下圖的狀況是須要刪除的:
刪除以後判斷新的前驅和後繼,接着刪除須要刪除的,直到前驅後繼合法爲止
效率:$O(nlogn)$(每一個點只會被插入和刪除最多一次,實際上用set,而後刪除時用迭代器遍歷可使得常數變得更小??)
這裏給出一個十分簡單的例題:
bzoj P2300==luogu P2521 [HAOI2011]防線修建
這個題目只要求維護上凸殼以及上凸殼的總長,並且由於題目的限制使得這個題的細節比模板簡單
用這個題的代碼代替模板的緣由是[數據刪除]
......好吧,是筆者set用得太爛了,給出的代碼實在很差意思當模板
#include<set> #include<cmath> #include<cstdio> #include<cstring> #include<algorithm> #define db double using namespace std; const db eps=1e-10; int dcmp(db x){ if(fabs(x)<eps)return 0; return x>0?1:-1; } struct Poi{ db x,y; Poi (db x=0,db y=0):x(x),y(y){ }; }; typedef Poi Vec; Vec operator + (Vec a,Vec b){ return Vec(a.x+b.x,a.y+b.y); } Vec operator - (Vec a,Vec b){ return Vec(a.x-b.x,a.y-b.y); } Vec operator * (Vec a,db t){ return Vec(a.x*t,a.y*t); } Vec operator * (db t,Vec a){ return Vec(a.x*t,a.y*t); } Vec operator / (Vec a,db t){ return Vec(a.x/t,a.y/t); } bool operator < (Vec a,Vec b){ return dcmp(a.x-b.x)<0||(!dcmp(a.x-b.x)&&dcmp(a.y-b.y)<0); }; bool operator > (Vec a,Vec b){ return dcmp(a.x-b.x)>0||(!dcmp(a.x-b.x)&&dcmp(a.y-b.y)>0); }; bool operator == (Vec a,Vec b){ return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y); }; set <Poi > C_Hull_set; set <int >aaa; int n,x,y,m,q,top; db an_stk[200010]; Vec city[100010],c0,c1,c2; bool lim[100010]; int ask[200010][2]; db Dot(Vec ,Vec ); db Cross(Vec ,Vec ); db Len(Vec ); db Area(Poi ,Poi ,Poi ); int drc(Poi ,Poi ,Poi ); db C_ins(Poi ,db ); int main() { int i,j,k; Poi now_l,now_r; scanf("%d%d%d",&n,&x,&y); c0=Poi(0,0),c1=Poi(x,y),c2=Poi(n,0); an_stk[0]=Len(c1-c0)+Len(c2-c1); C_Hull_set.insert(c0); C_Hull_set.insert(c1); C_Hull_set.insert(c2); scanf("%d",&m); for(i=1;i<=m;i++) scanf("%lf%lf",&city[i].x,&city[i].y); scanf("%d",&q); for(i=1;i<=q;i++){ scanf("%d",&ask[i][0]); if(ask[i][0]==1) scanf("%d",&ask[i][1]),lim[ask[i][1]]=true; } for(i=1;i<=m;i++) if(!lim[i]) an_stk[0]=C_ins(city[i],an_stk[0]); an_stk[++top]=an_stk[0]; for(i=q;i>0;i--){ if(ask[i][0]==2) an_stk[++top]=an_stk[top-1]; else an_stk[top]=C_ins(city[ask[i][1]],an_stk[top]); } for(i=top-1;i>0;i--) printf("%.2lf\n",an_stk[i]); return 0; } db Dot(Vec a,Vec b){ return a.x*b.x+a.y*b.y; } db Cross(Vec a,Vec b){ return a.x*b.y-b.x*a.y; } db Len(Vec v){ return sqrt(Dot(v,v)); } db Area(Poi a,Poi b,Poi c){ return Cross(b-a,c-a); } int drc(Poi p0,Poi p1,Poi p2){ return dcmp(Cross(p1-p0,p2-p0)); } db C_ins(Poi p,db las_ans){ set <Poi >:: iterator iter1; set <Poi >:: iterator iter2; set <Poi >:: iterator iter3; iter2=C_Hull_set.upper_bound(p); iter1=iter2,iter1--; if(Cross(p-*iter1,*iter2-*iter1)>=0)return las_ans; las_ans-=Len(*iter2-*iter1); iter2=iter1,iter2--; for(;!(*iter1==c0)&&drc(*iter2,*iter1,p)>=0;iter1--,iter2--,C_Hull_set.erase(iter3)){ las_ans-=Len(*iter1-*iter2); iter3=iter1; } iter1=C_Hull_set.upper_bound(p); iter2=iter1,iter2++; for(;!(*iter1==c2)&&drc(p,*iter1,*iter2)>=0;iter1++,iter2++,C_Hull_set.erase(iter3)){ las_ans-=Len(*iter1-*iter2); iter3=iter1; } iter2=C_Hull_set.upper_bound(p); iter1=iter2,iter1--; las_ans+=Len(*iter2-p)+Len(*iter1-p); C_Hull_set.insert(p); return las_ans; }
給直線加一個方向,他便成爲了有向直線,
有向直線能夠區分左右——逆時針轉動方向向量,率先掃過的一側爲左,另外一側爲右,
規定有向直線左側平面爲該有向直線肯定的半平面,
半平面的交集爲半平面交,
因爲咱們對直線的定義方式是線上一點與方向向量,
因此有向直線的定義只要沿用直線的定義,並確保更規範地使用方向向量便可
struct Line{ Point P; Vector v; Line (Point P,Vector v):P(P),v(v){ }; };//定義有向直線
bool operator == (Line l1,Line l2){ return !dcmp(atan2(l1.v.y,l1.v.x)-atan2(l2.v.y,l2.v.x)); }//定義直線比較爲極角比較 bool operator < (Line l1,Line l2){ return dcmp(atan2(l1.v.y,l1.v.x)-atan2(l2.v.y,l2.v.x))<0; }//定義直線比較爲極角比較 bool operator > (Line l1,Line l2){ return dcmp(atan2(l1.v.y,l1.v.x)-atan2(l2.v.y,l2.v.x))>0; }//定義直線比較爲極角比較
半平面是平面中的一個點集知足全部點都在某有向直線的左側,顯然是個凸集
因此半平面交也是個凸集
求半平面交的算法在咱們解決斜率優化時已經有所涉獵,
然而因爲斜率優化中自有的限制因此這裏給出的規範算法比以前的有更多的細節
這個算法與求凸包的方法相似:
排序,而後逐一掃描全部直線試圖將其加入表明半平面交邊界的隊列,在加入以前先肯定隊尾是否有直線會被徹底取代
方法是取出隊尾直線$l_i$和隊尾倒數第二條直線$l_{i-1}$判斷$l_i$與新加入的直線L的交點若他在$l_{i-1}$的右側則$l_i$能夠被剔除
重複此操做
須要注意的是因爲極角序是環形的,因此新加入的半平面可能繞了一圈把隊首的半平面剔除,因此還要額外判斷隊首是否須要被剔除
習題:
bzoj P1007==luogu P3194 因爲極角範圍使得這個題比模板簡單
圓還須要介紹嗎?
用圓心C和半徑r能夠肯定一個圓
引入參數方程能夠用一個弧度來肯定圓上一個點
圓的參數方程:
$$P=C+(r~cos\theta,r~sin\theta)$$
struct Circle{ Point c; double r; Circle (Point c,double r):c(c),r(r){ } Point poi(double cita){ return c+Point(cos(cita)*r,sin(cita)*r); } };//定義圓,引入圓的參數方程
(To be continue......)