計算幾何細節梳理&模板

點擊%XZY巨佬html

向量的板子

#include<bits/stdc++.h>
#define I inline
using namespace std;
typedef double DB;
struct Vec{
    DB x,y;
    I Vec(){x=y=0;}
    I Vec(DB a){x=a;y=0;}
    I Vec(DB a,DB b){x=a;y=b;}
    I friend istream&operator>>(istream&cin,Vec&a){return cin>>a.x>>a.y;}
    I friend ostream&operator<<(ostream&cou,Vec a){return cou<<a.x<<' '<<a.y;}
    I Vec operator-(){return Vec(-x,-y);}
    I Vec operator+(Vec a){return Vec(x+a.x,y+a.y);}
    I Vec operator-(Vec a){return Vec(x-a.x,y-a.y);}
    I Vec operator*(DB a){return Vec(x*a,y*a);}
    I Vec operator/(DB a){return Vec(x/a,y/a);}
    I friend Vec operator*(DB a,Vec b){b.x*=a,b.y*=a;return b;}
    I friend Vec operator/(DB a,Vec b){b.x/=a,b.y/=a;return b;}
    I Vec&operator+=(Vec a){x+=a.x,y+=a.y;return*this;}
    I Vec&operator-=(Vec a){x-=a.x,y-=a.y;return*this;}
    I Vec&operator*=(DB a){x*=a,y*=a;return*this;}
    I Vec&operator/=(DB a){x/=a,y/=a;return*this;}
    I DB operator*(Vec a){return x*a.x+y*a.y;}
    I DB operator^(Vec a){return x*a.y-y*a.x;}
    I friend bool operator<(Vec a,Vec b){
        return (a^b)>0||((a^b)==0&&Len(a)<Len(b));
    }
    I friend DB Len(Vec a){
        return sqrt(a.x*a.x+a.y*a.y);
    }
    I friend Vec Turn(Vec a,DB r){
        const DB c=cos(r),s=sin(r);
        return Vec(a.x*c-a.y*s,a.y*c+a.x*s);
    }
};

初階

向量運算

點積:\(x_1x_2+y_1y_2\),是一個向量在另外一個向量上的投影
叉積:\(x_1y_2-x_2y_1\),是兩個向量造成的平行四邊形的有向面積
用途很廣,搬一張圖ios

旋轉公式

\((x,y)\)\(r\)弧度\(\rightarrow(x\cos r-y\sin r,y\cos r+x\sin r)\)c++

點到直線距離

算一個叉積,除以底邊長就獲得了高,也就是點到直線距離。算法

DB Dis(Vec&a,Vec&b1,Vec&b2){
    return fabs((b1-a)^(b2-a))/Len(b1-b2);
}

直線交點&線段交點

點斜式要各類特判,不是理想的實現方法
考慮這樣的轉化方式:叉積->面積比->高的比->長度比->座標
以下圖,兩直線相交於\(P\)\(F\)\(G\)是垂足,有\(\triangle A_1PF\sim\triangle CB_1G\),能夠推出
\[\frac{\vec b×\vec c}{\vec b×\vec a}=\frac{S_{B_1B_2EA_1}}{S_{B_1B_2DC}}=\frac{|A_1F|}{|CG|}=\frac{|A_1P|}{|\vec a|}\]
\(P=A_1+\frac{|A_1P|}{|\vec a|}\vec{a}\),而後就作出來了。
蒟蒻到如今好像仍是背不了代碼,只能背下面這張圖了TAT
數組

如何判斷線段相交呢?先把直線交點求出來,再點積判一下位置關係不就好啦app

Vec LineCross(Vec a1,Vec a2,Vec b1,Vec b2){
    a2-=a1;b2-=b1;
    return b2^a2?a1+(b2^(b1-a1))/(b2^a2)*a2:Vec(NAN,NAN);
}
Vec SegCross(Vec&a1,Vec&a2,Vec&b1,Vec&b2){
    Vec c=LineCross(a1,a2,b1,b2);
    return (a1-c)*(a2-c)>0||(b1-c)*(b2-c)>0?Vec(NAN,NAN):c;
}

判斷點是否在多邊形內

過點任做一條射線(固然平行座標軸的最方便啦),與多邊形交奇數次則在多邊形內。dom

凸包

極角排序

能夠用atan2,也能夠直接用叉積判斷。
把最下面的點的座標找出來,把全部點總體平移,這時候全部點都不在\(x\)軸下面了,在\(180°\)內叉積能夠起到比較方向的做用。
代碼見凸包部分。ide

求凸包

單調棧維護。如圖,\((P_i-P_{t-1})×(P_t-P_{t-1})\ge0\),則\(P_t\)不在凸包上,須要彈掉。

洛谷P2742 【模板】二維凸包 / [USACO5.1]圈奶牛Fencing the Cows函數

I bool Polar(Vec&a,Vec&b){return (a^b)>0||((a^b)==0&&Len(a)<Len(b));}
//即Vec模板裏重載的小於號。共線向量把短的放前面,求凸包的時候方便彈掉
I int Convex(Vec*a,Vec*e){
    R n=e-a,k=0;
    for(R i=1;i<n;++i)
        if(a[i].y<a[k].y||(a[i].y==a[k].y&&a[i].x<a[k].x))k=i;
    swap(a[0],a[k]);
    const Vec tmp=a[0];
    for(R i=0;i<n;++i)a[i]-=tmp;
    sort(a+1,a+n,Polar);
    R*st=new int[n],p=0;
    for(R i=1;i<n;st[++p]=i++)
        while(p&&((a[i]-a[st[p-1]])^(a[st[p]]-a[st[p-1]]))>=0)--p;
    for(R i=0;i<=p;++i)a[i]=a[st[i]]+tmp;
    return p+1;
}
int main(){
    R n;cin>>n;
    for(R i=1;i<=n;++i)cin>>a[i];
    n=Convex(a+1,a+n+1);
    DB ans=Len(a[1]-a[n]);
    for(R i=1;i<n;++i)ans+=Len(a[i+1]-a[i]);
    printf("%.2lf\n",ans);
    return 0;
}

判斷點是否在凸包內

須要令凸包最下面的點爲\((0,0)\),不要還原。
把點的座標也作相應的變換以後,丟進數組裏lower_bound找到與它極角最接近的兩個點造成的線段,叉積判便可。ui

bool Inside(Vec*a,Vec*e,Vec v){
    R i=lower_bound(a,e,v)-a-1;
    return ((a[i+1]-a[i])^(v-a[i]))>=0;
}

多邊形面積

對於凸多邊形的理解:將多邊形劃分紅若干個三角形,每一個三角形的面積是兩個向量叉積的\(\frac12\)
update:其實任意多邊形面積均可以用三角剖分法求,詳見洛穀日報 計算幾何初步 by wjyyy

DB PolygonArea(Vec*a,Vec*e){
    R n=e-a;DB s=0;
    for(R i=2;i<n;++i)s+=(a[i-1]-a[0])^(a[i]-a[0]);
    return s/2;
}

進階算法

凸包

旋轉卡殼

利用凸包的單調性,在枚舉點的時候更新決策點(對踵點)並更新答案。
洛谷P1452 Beauty Contest

int main(){
    R n;cin>>n;
    for(R i=1;i<=n;++i)cin>>a[i];
    n=Convex(a+1,a+n+1);
    DB ans=0;
    a[++n]=a[1];//最後加一個點
    for(R i=1,p=1;i<n;++i){
        while(Dis(a[p],a[i],a[i+1])<Dis(a[p+1],a[i],a[i+1]))p=p%n+1;
        ans=max(ans,max(Len2(a[p]-a[i]),Len2(a[p]-a[i+1])));
    }
    printf("%.0lf\n",ans);
    return 0;
}

半平面交

XZY巨佬提到了一種\(O(n^2)\)的算法。其實,若是不是在線插入的話,能夠作到\(O(n\log n)\)
咱們用有向直線(一個點和一個方向向量)表示半平面,如下默認半平面在有向直線的左側。
對有向直線按方向向量的極角排序,維護一個雙端隊列,存儲當前構成半平面的直線以及相鄰兩直線的交點。
每次加入一條有向直線,若是隊首/隊尾的交點在直線右側(用叉積判)則彈掉隊首/隊尾的直線。
爲何這樣是對的呢?由於加入直線的單調性,因此要被彈出的直線必定在隊首或隊尾。感興趣的話能夠本身手畫一些例子來理解。
須要注意的細節:

  1. 加入直線時,先彈隊尾,再彈隊首。
  2. 最後還要檢查隊尾交點是否在隊首直線的右側,若是是也要彈掉。
  3. 特判平行直線,在右側的要彈掉。
  4. 若是題目給出的半平面不必定有限制邊界,則應該手動加入一個INF邊界。

另有洛谷P3222 [HNOI2012]射箭 的題解
模板題:洛谷P4196 [CQOI2006]凸多邊形

struct Line{
    Vec p,v;DB ang;
    I Line(){}
    I Line(Vec a,Vec b){p=a,v=b-a,ang=atan2(v.y,v.x);}
    I bool operator<(Line&a){return ang<a.ang;}
    I bool Right(Vec&a){return (v^(a-p))<-EPS;}
    I friend Vec Cross(Line&a,Line&b){return a.p+(b.v^(b.p-a.p))/(b.v^a.v)*a.v;}
}a[N],q[N];
DB HalfPlane(Line*a,Line*e){
    R n=e-a,h=0,t=0;
    sort(a,e);q[0]=a[0];
    for(R i=1;i<n;++i){
        while(h<t&&a[i].Right(k[t-1]))--t;
        while(h<t&&a[i].Right(k[h]))++h;
        if(q[t].ang!=a[i].ang)q[++t]=a[i];
        else if(a[i].Right(q[t].p))q[t]=a[i];
        if(h<t)k[t-1]=Cross(q[t-1],q[t]);
    }
    while(h<t&&q[h].Right(k[t-1]))--t;
    k[t]=Cross(q[t],q[h]);
    return PolygonArea(k+h,k+t+1);
}
int main(){
    R n,m,t=0;
    cin>>n;
    for(R i=1;i<=n;++i){
        Vec fst,lst,now;
        cin>>m>>fst;lst=fst;
        for(R j=2;j<=m;++j)
            cin>>now,a[++t]=Line(lst,now),lst=now;
        a[++t]=Line(lst,fst);
    }
    printf("%.3lf\n",HalfPlane(a+1,a+t+1));
    return 0;
}

閔可夫斯基和

對於歐氏空間的兩個點集\(A,B\),其閔可夫斯基和爲點集\(C=\{a+b|a\in A,b\in B\}\)
其中加法就是向量加法。
接下來咱們只討論凸包的閔可夫斯基和。

一個四邊形,一個三角形,兩個凸包的閔可夫斯基和會長什麼樣子呢?
由於凸包的特性,咱們只須要取其中一個凸集全部最外層的點,將另外一個凸多邊形沿這個點向量移動,就能夠獲得閔可夫斯基和。

觀察一下它的特色:外面正好有\(7\)個點,\(7\)條邊,有沒有注意到\(7=4+3\)呢?
實際上,對於任意兩個凸包來講,它們的閔可夫斯基和與它們的每一條邊按極角序順次相連所得的圖形都是全等的,也是一個凸包。
那豈不是很好求?二路歸併便可。須要處理三點共線的狀況,了撇一點的話直接扔進Convex函數裏再求一遍凸包就行了。
洛谷P4557 [JSOI2018]戰爭

typedef long long DB;//不涉及小數運算,直接開longlong
int Convex(Vec*a,Vec*e,Vec&bs){
    R n=e-a,k=0,p=0;
    for(R i=1;i<n;++i)
        if(a[i].y<a[k].y||(a[i].y==a[k].y&&a[i].x<a[k].x))k=i;
    swap(a[0],a[k]);bs+=a[0];
    for(R i=n-1;~i;--i)a[i]-=a[0];
    sort(a,e);
    for(R i=1;i<n;st[++p]=i++)
        while(p&&((a[i]-a[st[p-1]])^(a[st[p]]-a[st[p-1]]))>=0)--p;
    for(R i=0;i<=p;++i)a[i]=a[st[i]];//作出來沒還原
    return a[p+1]=0,p+1;
}
void Minkowski(R n,R m){
    for(R i=n;i;--i)a[i]-=a[i-1];
    for(R j=m;j;--j)b[j]-=b[j-1];
    for(R i=1,j=1,k=0;i<=n||j<=m;++k)
        c[k]=c[k-1]+(i<=n&&(j>m||a[i]<b[j])?a[i++]:b[j++]);
}
int main(){
    ios::sync_with_stdio(0);
    R n=in(),m=in(),q=in();
    Vec bs=0,v;
    for(R i=0;i<n;++i)cin>>a[i];
    for(R i=0;i<m;++i)cin>>b[i],b[i]=-b[i];
    n=Convex(a,a+n,bs);
    m=Convex(b,b+m,bs);
    Minkowski(n,m);
    n=Convex(c,c+n+m+1,v);
    while(q--)
        cin>>v,cout<<Inside(c,c+n,v-bs)<<'\n';
    return 0;
}

其它

最小圓覆蓋

把點集random_shuffle
第一層枚舉一個點,若是不在當前圓中,令其爲圓心,半徑設爲\(0\)
第二層枚舉另外一個點,若是不在當前圓中,將圓更新爲以當前兩個點的連線爲直徑的圓。
第三層枚舉另外一個點,若是不在當前圓中,將圓更新爲當前三個點的外接圓。
指望複雜度\(O(n)\)
洛谷P1742 最小圓覆蓋

typedef long double DB;//此題略卡精度
I void PerpLine(Vec a1,Vec a2,Vec&b1,Vec&b2){//中垂線
    b1=(a1+a2)/2;b2=b1+Turn(a1-a2);
}
I void CircleCoverage(Vec*a,Vec*e,Vec&O,DB&r){
    R n=e-a;
    for(R i=0;i<n;++i){
        if(Len(a[i]-O)<=r)continue;
        O=a[i],r=0;
        for(R j=0;j<i;++j){
            if(Len(a[j]-O)<=r)continue;
            O=(a[i]+a[j])/2,r=Len(a[i]-a[j])/2;
            Vec a2=O+Turn(a[i]-a[j]),b1,b2;
            for(R k=0;k<j;++k){
                if(Len(a[k]-O)<=r)continue;
                PerpLine(a[i],a[k],b1,b2);
                O=LineCross(O,a2,b1,b2),r=Len(a[i]-O);
            }
        }
    }
}
int main(){
    R n;cin>>n;
    for(R i=1;i<=n;++i)cin>>a[i];
    srand(20020307);
    random_shuffle(a+1,a+n+1);
    Vec O;DB r=0;
    CircleCoverage(a+1,a+n+1,O,r);
    cout<<fixed<<setprecision(10)<<r<<endl<<O<<endl;
    return 0;
}

自適應辛普森積分

求函數的定積分,應用於算不規則圖形的面積。採用二次函數來擬合。
\[\int_l^r f(x)dx\approx\int_l^r(ax^2+bx+c)dx\\ =\frac a3(r-l)^3+\frac b2(r-l)^2+c(r-l)\\ =\frac{r-l}6[2a(r^2-rl+l^2)+3b(r-l)+6c]\\ =\frac{r-l}6[f(l)+f(r)+4f(\frac{l+r}2)]\]
不斷分割區間並迭代,直到偏差控制在EPS之內。
由於偏差和區間長度有關,因此能夠將EPS也進行迭代,每次除以\(2\)
洛谷P4526 【模板】自適應辛普森法2

#include<bits/stdc++.h>
#define DB double
using namespace std;
DB a,EPS=1e-6;
inline DB f(DB x){
    return pow(x,a/x-x);
}
inline DB Calc(DB l,DB r){
    return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}
DB Simpson(DB l,DB r,DB ans,DB EPS){
    DB m=(l+r)/2,al=Calc(l,m),ar=Calc(m,r);
    return fabs(al+ar-ans)<EPS?ans:Simpson(l,m,al,EPS/2)+Simpson(m,r,ar,EPS/2);
}
int main(){
    cin>>a;
    if(a<0)puts("orz");
    else printf("%.5lf\n",Simpson(EPS,20,Calc(0,20),EPS));
    return 0;
}
相關文章
相關標籤/搜索