計算幾何相關
向量表示法
這里最適合的就是用一個二維點對\((a,b)\)來表示了。
點積
\({a.x*b.x+a.y*b.y}\)
在向量的含義下:\(\vec{a}·\vec{b}=|\vec{a}||\vec{b}|cos<\vec{a},\vec b>\)
叉積
\({a.x*b.y-a.y*b.x}\)
這個東西很有用,首先這個東西的絕對值就是兩個向量構成的三角形的面積的二倍。
證明的話只需要把圖畫出來,然后過向量端點的四條平行於坐標軸的直線將包含了整個向量的矩形分割成四個部分,然后割補法算一算就好了。
發現上面寫面積的時候帶了絕對值,當\(b\)在\(a\)的逆時針方向時,即\(b\)在\(a\)的左側時,兩者的叉積為正,否則叉積為負,所以通過叉積也可以很方便的計算向量之間的位置關系。
轉角
把向量轉動一個角度\(\alpha\)。
顯然向量如果把它的模長變為\(1\),就可以在單位圓上對應一個點,假設其模長為\(R\),那么一個向量可以寫為\((R\cos\beta,R\sin\beta)\),那么轉角只需要改變\(\beta\)即可。
再用三角函數的和角公式可以很容易拆分。
極角
其實就是如何計算上面那個\(\beta\),直接\(atan2(a.y,a.x)\)就好了。
兩條直線的交點
接下來所有的直線基本都用線段代替。
而線段基本都用一個點+一個方向向量表示。
假設\(a1,b1\)是點,\(a2,b2\)是方向向量。
\(b2\times (b1-a1))\)與\(b2\times a2\)的比值就只與交點到直線\(b\)與\(a1+a2\)到直線\(b\)的距離相關。
那么通過這個比例可以很容易算出交點的確切位置。
線段之間是否有交
似乎可以算下直線交點,判斷是否在線段上。
假設\(a_1,a_2,b_1,b_2\)是兩條線段的端點。
如果\((b_1-a_1)\times (b1-a2)\)與\((b_2-a_1)\times (b_2-a_2)\)異號,
並且\((a_1-b_1)\times (a_1-b_2)\)與\((a_2-b_1)\times (a_2-b_2)\)異號,
那么兩條線段有交。
上面的一部分代碼
struct Point{double x,y,ang;};
bool operator<(Point a,Point b){return (a.ang!=b.ang)?a.ang<b.ang:a.x<b.x;}
Point operator+(Point a,Point b){return (Point){a.x+b.x,a.y+b.y};}
Point operator-(Point a,Point b){return (Point){a.x-b.x,a.y-b.y};}
Point operator*(Point a,double b){return (Point){a.x*b,a.y*b};}
Point operator/(Point a,double b){return (Point){a.x/b,a.y/b};}
double operator*(Point a,Point b){return a.x*b.x+a.y*b.y;}
double Cross(Point a,Point b){return a.x*b.y-a.y*b.x;}
double Len(Point a){return sqrt(a.x*a.x+a.y*a.y);}
double Dis(Point a,Point b){return Len(a-b);}
Point Rotate(Point p,double a){double c=cos(a),s=sin(a);return (Point){p.x*c-p.y*s,p.x*s+p.y*c};}
struct Line{Point a,v;};
Point Intersection(Line a,Line b)
{
Point c=b.a-a.a;
double t=Cross(b.v,c)/Cross(b.v,a.v);
return a.a+a.v*t;
}
Graham掃描法
找到最左下角的點,其他點以這個點為原點進行極角排序。
按照次序依次考慮每個點加進凸包后是否能夠刪去之前的一些點。
通過叉積可以判斷當前凸包中的是否有可以被刪去的點。
時間復雜度\(O(nlogn)\),瓶頸在於極角排序。
叉積判斷那里判斷方法很多,並不唯一。
求凸包前可以給坐標加上\(\mbox{eps}\),防止一些毒瘤出題人造了一堆奇怪的數據導致凸包跑出來有問題。(比如說共線之類的)
Point S[MAX];int top;
void Graham(Point *p,int n)
{
int pos=1;
for(int i=2;i<=n;++i)
if(p[i].x<p[pos].x||(p[i].x==p[pos].x&&p[i].y<p[pos].y))
pos=i;
swap(p[1],p[pos]);
for(int i=2;i<=n;++i)p[i].ang=atan2(p[i].y-p[1].y,p[i].x-p[1].x);
sort(&p[2],&p[n+1]);S[++top]=p[1];S[++top]=p[2];
for(int i=3;i<=n;++i)
{
while(top>2&&Cross(p[i]-S[top],p[i]-S[top-1])>=0)--top;
S[++top]=p[i];
}
}
凸包面積
算出凸包之后選定一條邊,將凸包分割為若干三角形,叉積計算面積。
注意一下叉積的正負性,不要亂來。
判斷點是否在多邊形內
向一側做一條射線,計算與多邊形邊界的交的次數,如果為奇數則在內部,否則在外部。
注意射線過頂點的情況,最好給斜率加上一定擾動,防止這類情況的出現。
判斷點是否在凸包內
一種傻逼做法就是算這個點和凸包的每條邊形成的三角形面積和,看看是否等於凸包面積。
機智點的話,找到凸包上極角在這個點兩側的兩個點,叉積隨便判斷一下就好了。
旋轉卡殼(怎么讀我不負責)
比如說用來求解凸包的直徑,即凸包上最遠的兩點之間的距離。
通過凸包的單調性來實現。
比如求解凸包的直徑:
double Diameter(Point *S,int n)
{
double d=0;S[n+1]=S[1];if(n==2)return Dis(S[1],S[2]);
for(int i=1,j=3;i<=n;++i)
{
while(Cross(S[j]-S[i],S[j]-S[i+1])<Cross(S[j+1]-S[i],S[j+1]-S[i+1]))
j=(j==n)?1:j+1;
d=max(d,max(Dis(S[j],S[i]),Dis(S[j],S[i+1])));
}
return d;
}
所有要求的東西都會隨着你的逆時針旋轉而旋轉,所以直接利用單調性暴力向前即可。
這樣子的復雜度是線性的。
再比如說求解最小矩形覆蓋
void ScanLine(int n)
{
S[n+1]=S[1];S[0]=S[n];
for(int i=1,j1=3,j2=3,j3=n;i<=n;++i)
{
if(i==1)while((S[i]-S[i+1])*(S[j3-1]-S[j3])>0)j3=(j3==1)?n:j3-1;
while(Cross(S[j1]-S[i],S[j1]-S[i+1])<=Cross(S[j1+1]-S[i],S[j1+1]-S[i+1]))j1=(j1==n)?1:j1+1;
while((S[i+1]-S[i])*(S[j2+1]-S[j2])>0)j2=(j2==n)?1:j2+1;
while((S[i+1]-S[i])*(S[j3+1]-S[j3])<0)j3=(j3==n)?1:j3+1;
Line l0=(Line){S[i],S[i+1]-S[i]};
Line l1=(Line){S[j1],S[i]-S[i+1]};
Line l2=(Line){S[j2],Rotate(S[i+1]-S[i],Pi/2)};
Line l3=(Line){S[j3],Rotate(S[i]-S[i+1],Pi/2)};
tmp[1]=Intersection(l0,l2);tmp[2]=Intersection(l2,l1);
tmp[3]=Intersection(l1,l3);tmp[4]=Intersection(l3,l0);
double area=Dis(tmp[1],tmp[2])*Dis(tmp[2],tmp[3]);
if(area<ans)
ans=area,Ans[1]=tmp[1],Ans[2]=tmp[2],Ans[3]=tmp[3],Ans[4]=tmp[4];
}
}
注意到\(i=1\)的時候先順時針找到了正確的位置,接下來才能逆時針向前繼續尋找。
最小圓覆蓋
一個暴力的做法。
枚舉第一個點,考慮當前圓是否包含了這個點,如果沒有,則把圓變成以這個點為圓心,半徑為\(0\)的圓。再枚舉第二個點,考慮圓是否包含了這個點,如果沒有,則把圓變成以這兩個點的中點為圓心,半徑為兩點距離一半的圓。再枚舉第三個點,節點是否在圓內,如果不在,則把圓直接變成這三個點的外接圓。
看起來是三方的,然而隨機打亂點之后期望\(O(n)\)
代碼戳這里
半平面交
看起來似乎是有兩種方法的。
第一種是在線的增量法。
先存下當前的半平面,初始時我們認為是一個極大的凸多邊形。
每次新加入一個向量(代表一條直線),讓其和前面的所有邊求交,只有這向量左側的才是合法的東西,然后把新的交點給加入進來。時間復雜度是\(O(n^2)\)的。這種方法支持隨時增量,可以查詢中間任意時刻的半平面交。
代碼如下,注意是否有交時的正負號,以及是否有等於的情況。
void pre()
{
top=0;
S[++top]=(Point){inf,inf};
S[++top]=(Point){-inf,inf};
S[++top]=(Point){-inf,-inf};
S[++top]=(Point){inf,-inf};
}
void Cut(Line a)
{
S[top+1]=S[1];int tot=0;
for(int i=1;i<=top;++i)
{
double v1=Cross(a.v,S[i]-a.a);
double v2=Cross(a.v,S[i+1]-a.a);
if(v1>=0)tmp[++tot]=S[i];
if(v1*v2<0)tmp[++tot]=Intersection(a,(Line){S[i],S[i+1]-S[i]});
}
top=tot;for(int i=1;i<=top;++i)S[i]=tmp[i];
}
還有一種方法是\(O(nlogn)\)的,然而我不太會,所以先咕了。
自適應辛普森積分
求積分。
把函數\(F(x)\)看成二次函數去近似比,每次二分當前位置,如果左側估計出來的值加上右側估計出來值與整個區間估計出來的值之間的差小於\(eps\),那么就當做這個區間的值,否則則遞歸處理。
代碼如下:
double Simpson(double l,double r){return (r-l)*(F(l)+F(r)+4*F((l+r)/2))/6;}
double Simpson(double l,double r,double ans)
{
double mid=(l+r)/2,L=Simpson(l,mid),R=Simpson(mid,r);
if(fabs(L+R-ans)<eps)return ans;
return Simpson(l,mid,L)+Simpson(mid,r,R);
}
閔可夫斯基和
給定兩個點集\(A,B\),顯然這些點都可以用向量表示,它們的閔可夫斯基和\(C=\{\vec{a}+\vec{b}|\vec{a}\in A,\vec{b}\in B\}\)。
現在給定了兩個凸包,求他們的閔可夫斯基和。
畫圖之后不難發現他們的閔可夫斯基和形成的凸包其實就是這兩個凸包的所有邊按照極角排序的順序重新拼接在一起的。所以求出兩個凸包后直接按照逆時針順序歸並排序就可以得到閔可夫斯基和。
更加詳細的內容參考\(\mbox{lalaxu}\)博客
代碼:
void Minkowski(Point *p1,Point *p2,int n,int m)
{
p1[n+1]=p1[1];p2[m+1]=p2[1];
for(int i=1;i<=n;++i)t1[i]=p1[i+1]-p1[i];
for(int i=1;i<=m;++i)t2[i]=p2[i+1]-p2[i];
S[top=1]=p1[1]+p2[1];
int i=1,j=1;
while(i<=n&&j<=m)
{
++top;S[top]=S[top-1];
if(Cross(t1[i],t2[j])>=0)S[top]=S[top]+t1[i++];
else S[top]=S[top]+t2[j++];
}
while(i<=n)++top,S[top]=S[top-1]+t1[i++];
while(j<=m)++top,S[top]=S[top-1]+t2[j++];
}
注意一下我這份代碼求出來之后要把最后一個點給丟掉,要不然左下角的那個點會被計算兩次。
其他的
似乎沒有啥東西了?
三維計算幾何就算了,懶得搞了。