計算幾何基礎

計算幾何基礎

Tags：高級算法

---

前言

Noip爆炸後差點退組。
不取模設小上界的人大概不配學計算幾何吧。
copy一些網上的博客好了。
計算幾何詳解：https://blog.csdn.net/clover_hxy/article/details/53966405
凸包詳解：http://www.cnblogs.com/aiguona/p/7232243.html
旋轉卡殼：http://www.javashuo.com/article/p-coeckejb-na.html
題目標籤：WFJ
聲明一下：一下代碼片斷現已統一格式：*表示叉積，&表示點積，^表示數乘。符合博主如今的碼風。固然對於實數並不存在運算符混淆問題。

向量線段相關

點積

\(a.x*b.x+a.y*b.y\)

叉積

\(a.x*b.y-a.y*b.x\)

借圖：左邊是點積，右邊是叉積

轉角公式

逆時針旋轉角度\(B\)
原：\((cosAr,sinAr)\)
現：\((cos(A+B)r,sin(A+B)r)=((cosAcosB-sinAsinB)r,(sinAcosB+cosAsinB)r)\)
即：\((x,y)->(xcosB-ysinB,xsinB+ycosB)\)

極角排序

兩種方法

\(atan2\)

這是一個給定三角形對邊和鄰邊算角的函數，在cmath庫中
對每一個向量求\(angle=atan2(a.y,a.x)\)後直接排序便可

舉例：\(atan2(1,1)=0.78,atan2(1,0)=1.57\)

叉積

叉積的正負能夠判斷左右關係，因此在最大跨角不超過180度的時候能夠很是好地實現，不然有可能隨機一個起點繞好幾圈

兩直線夾角

若兩個向量都是從原點出發，則能夠很方便地用\[\theta =atan2(x*y,x\&y)\]叉積是平行四邊形，底乘高；點積是底乘投影。
更特殊點，該夾角是有向夾角，從x向量旋轉到y，逆時針是正角、順時針是負角。

兩直線交點

面積比->線段比->定比分點

Node Crosss(Node a1,Node a2,Node b1,Node b2)
{
    Node a=a2-a1,b=b2-b1,c=b1-a1;
    if(fabs(b*a)<1e-12) return (Node){-1e9,-1e9};//平行
    double t=(b*c)/(b*a);
    return a1+(a^t);
}

給個圖形象一點（以前畫錯了orzFlashHu）

判斷兩線段是否相交

\(a1,a2,b1,b2\)爲兩條直線的兩個端點
若是知足\((b1-a1)×(b1-a2)\)和\((b2-a1)×(b2-a2)\)不一樣號，而且\((a1-b1)×(a1-b2)\)和\((a2-b1)×(a2-b2)\)不一樣號，則兩線段相交，叫跨立實驗。

多邊形相關

凸包

找到最下的點設爲原點，將剩下的點用叉積極角排序
用單調棧維護凸包，當\((A[i]-A[sta[top-1]])×(A[sta[top]]-A[sta[top-1]])\)爲非負時須要彈棧

int cmp1(const Node&A,const Node&B) {return A.y<B.y||(A.y==B.y&&A.x<B.x);}
int cmp2(const Node&A,const Node&B) {return A*B>0||(A*B==0&&A.dis()<B.dis());}
//這個表示若是兩向量共線，把短的放前面（方便被踢開）
void Tubao()
{
    sort(A+1,A+n+1,cmp1);//找到最底下的點
    for(int i=n;i>=1;i--) A[i]=A[i]-A[1];
    sort(A+2,A+n+1,cmp2);//極角排序
    for(int i=1;i<=n;sta[++top]=i,i++)
        while(top>=2&&(A[i]-A[sta[top-1]])*(A[sta[top]]-A[sta[top-1]])>=0) top--;
    n=top;for(int i=1;i<=top;i++) A[i]=A[sta[i]];
}

判斷點是否在多邊形內

從該點向右引射線，與多邊形交奇數次即在其內。
順時針+1，逆時針-1，注意與頂點重合的狀況
\(Practice：Codeforces375C\ Circling\ Round\ Treasures\)

判斷點是否在凸包內

先斷定和邊界的關係
而後找到與其極角相鄰的兩點，憑此判斷
須保證\(A[1]=(0,0)\)

ll in(Node a)
{
    if(a*A[1]>0||A[tot]*a>0) return 0;
    ll ps=lower_bound(A+1,A+tot+1,a,cmp2)-A-1;
    return (a-A[ps])*(A[ps%tot+1]-A[ps])<=0;
}

求多邊形面積

逆時針極角排序後直接用叉積求

double CalcS()
{
    double res=0;
    for(int i=2;i<node;i++) res+=(A[i]-A[1])*(A[i+1]-A[1]);
    return res/2;
}

三角形重心

三點各座標的平均值\((\frac{x1+x2+x3}{3},\frac{y1+y2+y3}{3})\)

算法

半平面交

用來求解一類線性規劃問題？

首先存一個極大的凸多邊形，考慮增長一條直線
那麼將其視爲一個向量，左邊的爲合法半平面
複雜度\(O(n^2)\)

void Cut(Node a,Node b)//增長一個向量
{
    A[tt+1]=A[1];c=0;
    for(int i=1;i<=tt;i++)//枚舉每一條邊
    {
        double v1=(A[i]-a)*(A[i]-b);
        double v2=(A[i+1]-a)*(A[i+1]-b);
        if(v1>=0) C[++c]=A[i];//若是A[i]在合法平面內則存入
        if(v1*v2<0) Add(A[i],A[i+1],a,b);//若是邊(A[i],A[i+1])穿過新增直線，把交點加入
    }
    tt=c;for(int i=1;i<=tt;i++) A[i]=C[i];
}

若是不是凸多邊形的時候，能夠採用水平可見直線的方法，先按照斜率排序後直接用棧維護下平面。複雜度爲排序複雜度，有必定侷限性，但在座標範圍較大采用上面方法會炸精度時很是好用

struct Line{int k,b;int id;}A[N];
int cmp1(const Line&A,const Line&B) {return A.k<B.k||(A.k==B.k&&A.b>B.b);}
double Node(int i,int j) {return 1.0*(A[j].b-A[i].b)/(A[i].k-A[j].k);}
void work()
{
    sort(A+1,A+n+1,cmp1);
    for(int i=1;i<=n;sta[++top]=i,i++)
        while(top>=2&&Node(i,sta[top])-Node(sta[top],sta[top-1])<=eps) top--;
}

旋轉卡殼

利用凸包上一些奇妙的單調性，求解

• 多邊形直徑
• 多邊形寬度
• 最小面積矩形覆蓋
• 等等

複雜度\(O(n)\)

for(int i=1,p=1;i<=n;i++)
{
    //這份代碼只卡了一個點，還能夠同時卡多個點
    while(H(A[i],A[i+1],A[p%n+1])>=H(A[i],A[i+1],A[p])) p=p%n+1;
    Ans=max(Ans,max((A[p]-A[i]).dis(),(A[p]-A[i+1]).dis()));
}

最小圓覆蓋

隨機增量構造，複雜度\(O(n)\)
感受能夠被扁平三角形卡掉，可是實際上並不會，因爲\(ijk\)的枚舉順序，三角形必定是在嘗試把三條邊都做爲直徑都失敗後、纔會有三角形的外接圓。（2019.2.23 byGXY）

pair<Node,Node> Line(int k,int i)//中垂線
{
    Node a=A[k]-A[i];a.rot();
    Node c=(A[k]+A[i])^0.5;
    return mp(c,a+c);
}
int main()
{
    random_shuffle(A+1,A+n+1);//先隨機化
    for(int i=1,j,k;i<=n;i++)
        if((A[i]-C).dis()>R)//某點不在圓內，更新爲半徑爲0的圓
            for(C=A[i],R=0,j=1;j<i;j++)
                if((A[j]-C).dis()>R)//某線段不在圓內，更新爲以該線段爲直線的圓
                    for(C=(A[i]+A[j])^0.5,R=(A[i]-C).dis(),k=1;k<j;k++)
                        if((A[k]-C).dis()>R)//某三角形不在圓內，更新爲該三角形的外切圓
                            C=Cross(Line(k,i),Line(k,j)),R=(A[i]-C).dis();//求兩條中垂線交點
}

自適應辛普森積分

用來求定積分。
當函數極其鬼畜、退不出原函數的時候，採用二次函數擬合，在精度達到要求以後便視爲相等
即\[\int_L^Rf(x)dx\approx\int_L^R(Ax^2+Bx+C)dx=\frac{(R-L)(f(L)+f(R)+4f(\frac{L+R}{2}))}{6}\]

double f(double x) {return (c*x+d)/(a*x+b);}
double calc(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,ans1=calc(l,mid),ans2=calc(mid,r);
    if(fabs(ans1+ans2-ans)<=eps) return ans;
    return simpson(l,mid,ans1)+simpson(mid,r,ans2);
}
int main()
{
    std::cin>>a>>b>>c>>d>>L>>R;
    printf("%.6f\n",simpson(L,R,calc(L,R)));
}

注意\(\frac{1}{x}\)的原函數是\(ln|x|\)！！

三維凸包

出門左拐：http://www.javashuo.com/article/p-pzldgser-kd.html

閔可夫斯基和

出門左拐第二家：http://www.javashuo.com/article/p-wkxofllf-kg.html

Pick定理、歐拉公式和圓的反演

出門左拐第三家：http://www.javashuo.com/article/p-rvbjiaoo-kh.html

Voronoi圖與Delaunay三角剖分

出門左拐第四家：http://www.javashuo.com/article/p-mwazquwp-kh.html

部分好題題解

[ZJOI2018] 保鏢