凸包

凸包

凸包(Convex Hull)是一個計算幾何(圖形學)中的概念。

在一個實數向量空間V中,對於給定集合X,所有包含X的凸集交集S被稱為X的凸包。X的凸包可以用X內所有點(X1,...Xn)的凸組合來構造.

在二維歐幾里得空間中,凸包可想像為一條剛好包著所有點的橡皮圈。

用不嚴謹的話來講,給定二維平面上的點集,凸包就是將最外層的點連線起來構成的凸多邊形,它能包含點集中所有的點。

基本介紹

  • 中文名:凸包
  • 外文名:Convex Hull
  • 定義:包含集合X的所有凸集交集
  • 內容:一個計算幾何(圖形學)中的概念
  • 條件:n維歐幾里得空間
定義,平面求法,代碼實例,

定義

⒈對於一個集合D,D中任意有限個點的凸組合的全體稱為D的凸包。
⒉對於一個集合D,所有包含D的凸集之交稱為D的凸包。
可以證明,上述兩種定義是等價的
概念
1  點集Q的凸包(convex hull)是指一個最小凸多邊形,滿足Q中的點或者在多邊形邊上或者在其內。右圖中由紅色線段表示的多邊形就是點集Q={p0,p1,...p12}的凸包。
示例圖(一)示例圖(一)
2  一組平面上的點,求一個包含所有點的最小的凸多邊形,這就是凸包問題了。這可以形象地想成這樣:在地上放置一些不可移動的木樁,用一根繩子把他們儘量緊地圈起來,並且為凸邊形,這就是凸包了。
數學定義:設S為歐幾里得空間
的任意子集。包含S的最小凸集稱為S的凸包,記作conv(S)。

平面求法

常見求法
凸包最常用的凸包算法是Graham掃描法和Jarvis步進法
Graham's Scan法
這個算法是由數學大師葛立恆(Graham)發明的,他曾經是美國數學學會(AMS)主席、AT&T首席科學家以及國際雜技師協會(IJA)主席。
問題
給定平面上的二維點集,求解其凸包。
過程
⒈ 在所有點中選取y坐標最小的一點H,當作基點。如果存在多個點的y坐標都為最小值,則選取x坐標最小的一點。坐標相同的點應排除。然後按照其它各點p和基點構成的向量<H,p>;與x軸的夾角進行排序,夾角由大至小進行順時針掃描,反之則進行逆時針掃描。實現中無需求得夾角,只需根據餘弦定理求出向量夾角的餘弦值即可。以下圖為例,基點為H,根據夾角由小至大排序後依次為H,K,C,D,L,F,G,E,I,B,A,J。下面進行逆時針掃描。
示例圖(二)示例圖(二)
⒉ 線段<H,K>;一定在凸包上,接著加入C。假設線段<K,C>;也在凸包上,因為就H,K,C三點而言,它們的凸包就是由此三點所組成。但是接下來加入D時會發現,線段<K,D>;才會在凸包上,所以將線段<K,C>;排除,C點不可能是凸包。
⒊ 即當加入一點時,必須考慮到前面的線段是否在凸包上。從基點開始,凸包上每條相臨的線段的旋轉方向應該一致,並與掃描的方向相反。如果發現新加的點使得新線段與上線段的旋轉方向發生變化,則可判定上一點必然不在凸包上。實現時可用向量叉積進行判斷,設新加入的點為
,上一點為
,再上一點為
。順時針掃描時,如果向量
的叉積為正(逆時針掃描判斷是否為負),則將上一點刪除。刪除過程需要回溯,將之前所有叉積符號相反的點都刪除,然後將新點加入凸包。
在上圖中,加入K點時,由於線段<H,C>要旋轉到<H,K>的角度,為順時針旋轉,所以C點不在凸包上,應該刪除,保留K點。接著加入D點,由於線段<K,D>要旋轉到<H,K>的角度,為逆時針旋轉,故D點保留。按照上述步驟進行掃描,直到點集中所有的點都遍歷完成,即得到凸包。
這個算法可以直接在原數據上進行運算,因此空間複雜度為O⑴。但如果將凸包的結果存儲到另一數組中,則可能在代碼級別進行最佳化。由於在掃描凸包前要進行排序,因此時間複雜度至少為快速排序的O(nlgn)。後面的掃描過程複雜度為O(n),因此整個算法的複雜度為O(nlgn)。
Jarvis步進法
對於一個有三個或以上點的點集Q,過程如下:
計算點集最右邊的點為凸包的頂點的起點,如上圖的P3點。
Do
For i = 0 To 總頂點數
計算有向向量P3->Pi
If 其餘頂點全部在有向向量P3->Pi的左側或右側,則Pi點為凸包的下一頂點
Pi點加入凸包列表
GoTo 1
End If
Next
Exit Do
1:
Loop
此過程執行後,點按極角自動順時針或逆時針排序,只需要按任意兩點的次序就可以了。而左側或右側的判斷可以用前述的矢量點積性質實現。
中心法
先構造一個中心點,然後將它與各點連線起來,按斜率遞增的方法,求出凸包上部;再按斜率遞減的方法,求出凸包下部。
水平法
從最左邊的點開始,按斜率遞增的方法,求出凸包上部;再按斜率遞減的方法,求出凸包下部。水平法較中心法減少了斜率無限大的可能,減少了代碼的複雜度
快包法
選擇最左、最右、最上、最下的點,它們必組成一個凸四邊形(或三角形)。這個四邊形內的點必定不在凸包上。然後將其餘的點按最接近的邊分成四部分,再進行快包法(QuickHull)。

代碼實例

C代碼一
#include<stdio.h>#include<math.h>#include<stdlib.h>typedefstruct{    doublex;    doubley;}POINT;POINTresult[102];//保存凸包上的點,相當於所說的棧SPOINTa[102];intn,top;doubleDistance(POINTp1,POINTp2)//兩點間的距離{    returnsqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));}doubleMultiply(POINTp1,POINTp2,POINTp3)//叉積{    return((p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x));}intCompare(constvoid*p1,constvoid*p2)//根據p0->p1的極值和p0->p2的極值進行比較,如果極值相同則用距離長度比較{    POINT*p3,*p4;    doublem;    p3=(POINT*)p1;    p4=(POINT*)p2;    m=Multiply(a[0],*p3,*p4);    if(m<0)return1;    elseif(m==0&&(Distance(a[0],*p3)<Distance(a[0],*p4)))        return1;    elsereturn-1;}//尋找凸包的過程,p0,p1,p2..的尋找過程在下面main中進行了voidTubao(){    inti;    result[0].x=a[0].x;    result[0].y=a[0].y;    result[1].x=a[1].x;    result[1].y=a[1].y;    result[2].x=a[2].x;    result[2].y=a[2].y;    top=2;    for(i=3;i<=n;i++)    {        while(Multiply(result[top-1],result[top],a[i])<=0&&top>2)            top--;        result[top+1].x=a[i].x;        result[top+1].y=a[i].y;        top++;    }}intmain(){    inti,p;    doublepx,py,len,temp;    while(scanf("%d",&n)!=EOF,n)    {        for(i=0;i<n;i++)            scanf("%lf%lf",&a[i].x,&a[i].y);        if(n==1)        {            printf("0.00/n");            continue;        }        elseif(n==2)        {            printf("%.2lf/n",Distance(a[0],a[1]));            continue;        }        //這裡的目的好像是找出y坐標最小的點,之後把他定義為P0        py=-1;        for(i=0;i<n;i++)        {            if(py==-1||a[i].y<py)            {                px=a[i].x;                py=a[i].y;                p=i;            }            elseif(a[i].y==py&&a[i].x<px)            {                px=a[i].x;                py=a[i].y;                p=i;            }        }        //swap(a[0],a[p])        temp=a[0].x;        a[0].x=a[p].x;        a[p].x=temp;        temp=a[0].y;        a[0].y=a[p].y;        a[p].y=temp;        //用叉乘來實現排序的比較        qsort(&a[1],n-1,sizeof(double)*2,Compare);        a[n].x=a[0].x;        a[n].y=a[0].y;        //調用tubao()        Tubao();        len=0.0;        for(i=0;i<top;i++)            len=len+Distance(result[i],result[i+1]);        printf("%.2lf/n",len);    }    return0;}
C代碼二
#include<iostream>//求點集合的凸包的gram算法。n是頂點個數,x,y是頂點坐標。#include<fstream>//order是按照頂點和左下腳的角度的排序後數組。#include<deque>//tu即是逆時針的凸包上的頂點。#include<math.h>usingnamespacestd;//使用條件:1。點可以任意給,可重複。//2。三個以及以上的點。ifstreamfin("input.txt");//3。已經考慮了邊上有點的情況。#defineNN1000#definepi3.1415827typedefstructCseg{    doublex,y,tg;}Cseg;intn;doublex[NN],y[NN];deque<Cseg>order;deque<int>tu;Csegseg1;deque<Cseg>::iteratorp1;deque<int>::iteratorp,q;voidin();voidgram();voidmakeorder(ints);doubledist(doublex1,doubleyy1,doublex2,doubleyy2);doublecross(doublex1,doubleyy1,doublex2,doubleyy2);voidout();intmain(){    in();    gram();    out();    return0;}voidout(){    for(inti=0;i<tu.size();i++)        cout<<order[tu].x<<""<<order[tu].y<<endl;    cout<<tu.size()<<"EdgesPolydgon"<<endl;    return;}voidin(){    fin>>n;    for(inti=0;i<n;i++)        fin>>x>>y;    return;}voidgram(){    intmm;    mm=0;    for(inti=1;i<n;i++)        if(y[mm]>y+1e-9)            mm=i;    elseif(fabs(y[mm]-y)<1e-9&&x[mm]>x+1e-9)        mm=i;    makeorder(mm);    seg1.x=x[mm];    seg1.y=y[mm];    tu.clear();    tu.push_back(0);    tu.push_back(1);    tu.push_back(2);    for(inti=3;i<order.size();i++)    {        p=tu.end();        seg1.x=order.x;        seg1.y=order.y;        p--;        q=p-1;        if(cross(order[*p].x-order[*q].x,order[*p].y-order[*q].y,order.x-order[*q].x,order.y-order[*q].y)>1e-9)            tu.push_back(i);        else        {            tu.pop_back();            i--;            continue;        }    }    return;}voidmakeorder(ints){    doubletg;    order.clear();    for(inti=0;i<n;i++)    {        if(i==s)            continue;        tg=atan2(y-y[s],x-x[s]);        seg1.x=x;        seg1.y=y;        seg1.tg=tg;        p1=order.begin();        while(p1!=order.end())        {            if(fabs(tg-p1->tg)<1e-9)            {                if(dist(x[s],y[s],x,y)>dist(x[s],y[s],p1->x,p1->y)+1e-9)                {                    p1->x=x;                    p1->y=y;                }                break;            }            elseif(tg<p1->tg)            {                order.insert(p1,seg1);                break;            }            p1++;        }        if(p1==order.end())            order.insert(p1,seg1);    }    seg1.x=x[s];    seg1.y=y[s];    order.insert(order.begin(),seg1);    return;}doublecross(doublex1,doubleyy1,doublex2,doubleyy2){    return(x1*yy2-x2*yy1);}doubledist(doublex1,doubleyy1,doublex2,doubleyy2){    returnpow((x1-x2)*(x1-x2)+(yy1-yy2)*(yy1-yy2),0.5);}
Pascal代碼三
const  pi=3.1415926575;  zero=1e-6;  maxn=1000;  maxnum=100000000;var  ans,temp:extended;  n,tot:longint;  x,y:array[0..maxn] of extended;  zz,num:array[0..maxn]of longint;procedure swap(var ii,jj:extended); var  t:extended;    begin      t:=ii;ii:=jj;jj:=t;    end;procedure init;var  i,j:longint;    begin      readln(n);       for i:=1 to n do readln(x[i],y[i]);   end;function ok(x,midx,y,midy:extended):longint;beginif  abs(x-midx)<=zero then  begin    if abs(midy-y)<=zero then exit(0);    if midy>y then exit(1)    else exit(2);  endelse   begin     if x<midx then exit(1)     else exit (2);   end;end;procedure qsort(head,tail:longint);var  i,j:longint;  midx,midy:extended;    begin      i:=head;      j:=tail;      midx:=x[(head+tail)div 2];      midy:=y[(head+tail)div 2];      repeat        while ok(x[i],midx,y[i],midy)=1 do inc(i);        while ok(x[j],midx,y[j],midy)=2 do dec(j);        if i<=j then            begin              swap(x[i],x[j]);              swap(y[i],y[j]);              inc(i);              dec(j);            end;     until i>j;      if i<tail then qsort(i,tail);      if j>head then qsort(head,j);end;function Plot(x1,y1,x2,y2:extended):extended;  begin     Plot:=x1*y2-x2*y1;  end;function check(first,last,new:longint):boolean;varax,ay,bx,by:extended;Pt:extended;begin  ax:=x[last]-x[first];ay:=y[last]-y[first];  bx:=x[new]-x[first];by:=y[new]-y[first];  if Plot(ax,ay,bx,by)<=0 then exit(true)   else exit(false);end;procedure Tbao;vari,j,tail:longint;begin tot:=0; zz[1]:=1;tail:=1; for i:=2 to n do    begin      while(zz[tail]<>1)and check(zz[tail-1],zz[tail],i) do dec(tail);      inc(tail);      zz[tail]:=i;    end;    inc(tot,tail-1);  for i:=1 to tail-1 do   num[i]:=zz[i];  zz[1]:=n;tail:=1;  for i:=n-1 downto 1 do  begin    while (zz[tail]<>n) and check(zz[tail-1],zz[tail],i) do dec(tail);    inc(tail);    zz[tail]:=i;  end;  for i:=1 to tail-1 do  num[tot+i]:=zz[i];  inc(tot,tail-1);end;function dist(a,b:longint):extended;   begin     dist:=sqrt((x[a]-x[b])*(x[a]-x[b])+(y[a]-y[b])*(y[a]-y[b]));   end;   procedure main; var    i,j:longint;  begin    qsort(1,n);    Tbao;    ans:=0;   for i:=1 to tot-1 do        ans:=ans+dist(num[i],num[i+1]);  ans:=ans+dist(num[tot],num[1]);  ans:=ans+temp*pi*2;  writeln(ans:0:1); end; begininit;main;end.

相關詞條

熱門詞條

聯絡我們