1. 程式人生 > >計算幾何基礎入門(1)

計算幾何基礎入門(1)

                                       平面最接近點對&&二維凸包&&最小覆蓋圓  模板

  首先要知道兩個基礎知識   叉積與基礎運算子過載   二維叉積可以用來判斷點與點的位置關係與面積  (三維叉積可算平面的法向量)

  這些模板都是二維的  我想想看能不能化為三維的題 空間最小點對 三維凸包 最小覆蓋球

  例題   P4894、 P3744、P2785

   

 

        平面最接近點對模板講解

 

 

 最容易想到的方法是n2暴力法  但一分析資料範圍就會發現肯定會t  於是我們有了下面一個思路  先將點按x排序再運用歸併把數對最終分為左右然後在合併   那麼最短距離就變為了左右最短距離與跨越分界線的最短距離這三種情況      難點就只在跨越分界線的了    關鍵點就是能夠點的距離小於d的點的數量是很少的  我們的任務就是隻找出並計算這些點  關鍵在程式碼內有說明

#include<bits/stdc++.h>
using namespace std;
const int maxn=1000001 ;
const int INF=2<<20;
int n,temp[maxn];
struct Point{double x,y;}S[maxn];
bool cmp(const Point &a,const Point&b ){if(a.x==b.x)return a.y<b.y; else return a.x<b.x;}
bool cmps(const int &a,const int &b){return S[a].y<S[b].y ;}
double min(double a,double b){return a<b?a:b;}
double dist(int i,int j) ///求距離
{
    double x=(S[i].x-S[j].x)*(S[i].x-S[j].x);
    double y=(S[i].y-S[j].y)*(S[i].y-S[j].y);
    return sqrt(x+y);
}

double merge(int left,int right) ///只判斷x小於d&&y也小於d的點對
{
    double d=INF;
    if(left==right)
        return d ;
    if(left+1==right)
        return dist(left,right);
    int mid=left+right>>1;
    double d1=merge(left,mid) ;///分
    double d2=merge(mid+1,right) ;
    d=min(d1,d2);
    int i,j,k=0;
    for(i=left;i<=right;i++)
        if(fabs(S[mid].x-S[i].x)<=d)///選出 點與分界線x距離小於d的點
            temp[k++]=i;
    sort(temp,temp+k,cmps);///按y排序
    for(i=0;i<k;i++)
        for(j=i+1;j<k&&S[temp[j]].y-S[temp[i]].y<d;j++)///選出點與點y距離小於d的點
        {
            double d3=dist(temp[i],temp[j]);///最終可證明每次只要判斷不超過6個點的距離  一個d*2d矩陣範圍內的點
            if(d>d3)d=d3;
        }
    return d;
}

int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%lf%lf",&S[i].x,&S[i].y);
    sort(S,S+n,cmp);///先按x排序 好劃分分界線
    printf("%.4lf\n",merge(0,n-1));
    return 0;
}

 

 二維凸包模板講解

 

 經分析只要我們選最外層的點組合成圖形就是最短圍欄  這就是二維凸包   這時我們採用的方法是Graham法  一個經典演算法   關鍵就在於如何判斷一個點是不是凸包的構成點   建議畫一個凸包並看程式碼分析   程式碼註釋寫了關鍵

#include<bits/stdc++.h>
using namespace std;
int n;
struct ben
{
    double x,y;
}p[10005],s[10005];
double check(ben a1,ben a2,ben b1,ben b2)///檢查叉積是否大於0,如果是a就逆時針轉到b
{
    return (a2.x-a1.x)*(b2.y-b1.y)-(b2.x-b1.x)*(a2.y-a1.y);
}
double d(ben p1,ben p2)///兩點間距離
{
    return sqrt((p2.y-p1.y)*(p2.y-p1.y)+(p2.x-p1.x)*(p2.x-p1.x));
}
bool cmp(ben p1,ben p2)///排序函式,這個函式別寫錯了,要不然功虧一簣
{
    double tmp=check(p[1],p1,p[1],p2);
    if(tmp>0)
        return 1;
    if(tmp==0&&d(p[0],p1)<d(p[0],p2))
        return 1;
    return 0;
}
int main()
{

    scanf("%d",&n);
    double mid;
    for(int i=1;i<=n;i++)
    {
        scanf("%lf%lf",&p[i].x,&p[i].y);
        if(i!=1&&p[i].y<p[1].y)///點的交換 p[1]為最低點
        {
            mid=p[1].y;p[1].y=p[i].y;p[i].y=mid;
            mid=p[1].x;p[1].x=p[i].x;p[i].x=mid;
        }
    }
    sort(p+2,p+1+n,cmp);///系統快排  利用叉積排序
    s[1]=p[1];///最低點一定在凸包裡
    int cnt=1;
    for(int i=2;i<=n;i++)///s陣列儲存凸包端點
    {                                                       ///怎麼判斷一個點會不會被踢走
                                                            /// 主要看點形成的邊是向右轉了還是向左轉了
        while(cnt>1&&check(s[cnt-1],s[cnt],s[cnt],p[i])<=0) ///判斷前面的會不會被踢走,如果被踢走那麼出棧
            cnt--;
        cnt++;
        s[cnt]=p[i];
    }
    s[cnt+1]=p[1];///最後一個點回到凸包起點
    double ans=0;
    for(int i=1;i<=cnt;i++)
        ans+=d(s[i],s[i+1]);///最小長度為邊長和
    printf("%.2lf\n",ans);
    return 0;
}

最小覆蓋圓模板講解

 

 這題的關鍵是  不斷向現在形成的圓加入新的點  如果這點不在此圓內時   那麼符合最小圓概念的話 這個新點它一定在符合的圓的邊界上   運用這個思路我們可以求出三個點    來求這個圓    由中垂線求圓心  中點直接/2  方向就是右轉90°   用向量來表示直線   

#include<bits/stdc++.h>
#define ll long long
using namespace std;

struct vec
{
    double x, y;
    vec (const double& x0 = 0, const double& y0 = 0) : x(x0), y(y0) {}///過載+-*/
    vec operator + (const vec& t) const {return vec(x+t.x, y+t.y);}
    vec operator - (const vec& t) const {return vec(x-t.x, y-t.y);}
    vec operator * (const double& t) const {return vec(x*t, y*t);}
    vec operator / (const double& t) const {return vec(x/t, y/t);}
    const double len2 () const {return x*x + y*y;}
    const double len () const {return sqrt(len2());}///向量模長
    vec norm() const {return *this/len();}///標準向量
    vec rotate_90_c () {return vec(y, -x);}///向右轉90度
};

double dot(const vec& a, const vec& b) {return a.x*b.x + a.y*b.y;}
double crs(const vec& a, const vec& b) {return a.x*b.y - a.y*b.x;}///叉積

vec lin_lin_int(const vec& p0, const vec& v0, const vec& p1, const vec& v1)
{
    double t = crs(p1-p0, v1) / crs(v0, v1);///運用叉積求交點
    return p0 + v0 * t;
}

vec circle(const vec& a, const vec& b, const vec& c)///中垂線求圓心  直線用向量表示
{
    return lin_lin_int((a+b)/2, (b-a).rotate_90_c(), (a+c)/2, (c-a).rotate_90_c());
}

int n;
vec pot[100005];

int main()
{
    scanf("%d", &n);
    for(int i=1; i<=n; i++) scanf("%lf%lf", &pot[i].x, &pot[i].y);
    random_shuffle(pot+1, pot+n+1);///隨機化  不隨機化會有大問題
    vec o;///o 為圓心  r2為半徑
    double r2 = 0;
    for(int i=1; i<=n; i++)
    {
        if((pot[i]-o).len2() > r2)///出現一個點不在現在圓的範圍內   立刻形成一個新圓
        {                         ///經分析 可知這個點一定在圓邊界上
            o = pot[i], r2 = 0;
            for(int j=1; j<i; j++)
            {
                if((pot[j]-o).len2() > r2)///類似  這時已確定兩個在邊界上的點了
                {
                    o = (pot[i]+pot[j])/2, r2 = (pot[j]-o).len2();
                    for(int k=1; k<j; k++)
                    {
                        if((pot[k]-o).len2() > r2)///三點確定圓
                        {
                            o = circle(pot[i], pot[j], pot[k]), r2 = (pot[k]-o).len2();
                        }
                    }
                }
            }
        }
    }
    printf("%.10lf\n%.10lf %.10lf\n", sqrt(r2), o.x, o.y);
    return 0;
}

 

看到最後面的彩蛋  : 輕音樂   Speak Softly,Love  &n