-
凸包
鎖定
凸包（Convex Hull）是一個計算幾何（圖形學）中的概念。
在二維歐幾里得空間中，凸包可想象為一條剛好包著所有點的橡皮圈。
用不嚴謹的話來講，給定二維平面上的點集，凸包就是將最外層的點連接起來構成的凸多邊形，它能包含點集中所有的點。
凸包定義
⒈對於一個集合D，D中任意有限個點的凸組合的全體稱為D的凸包。
⒉對於一個集合D，所有包含D的凸集之交稱為D的凸包。
可以證明，上述兩種定義是等價的
概念
圖1 示例圖
2 一組平面上的點，求一個包含所有點的最小的凸多邊形，這就是凸包問題了。這可以形象地想成這樣：在地上放置一些不可移動的木樁，用一根繩子把他們儘量緊地圈起來，並且為凸邊形，這就是凸包了。
凸包平面求法
常見求法
凸包最常用的凸包算法是Graham掃描法和Jarvis步進法
Graham's Scan法
這個算法是由數學大師葛立恆（Graham）發明的，他曾經是美國數學學會（AMS）主席、AT&T首席科學家以及國際雜技師協會（IJA）主席。
問題
給定平面上的二維點集，求解其凸包。
過程
⒈ 在所有點中選取y座標最小的一點H，當作基點。如果存在多個點的y座標都為最小值，則選取x座標最小的一點。座標相同的點應排除。然後按照其它各點p和基點構成的向量<H,p>；與x軸的夾角進行排序，夾角由大至小進行順時針掃描，反之則進行逆時針掃描。實現中無需求得夾角，只需根據餘弦定理求出向量夾角的餘弦值即可。以圖2為例，基點為H，根據夾角由小至大排序後依次為H，K，C，D，L，F，G，E，I，B，A，J。下面進行逆時針掃描。
圖2 示例圖
⒉ 線段<H,K>；一定在凸包上，接着加入C。假設線段<K,C>；也在凸包上，因為就H，K，C三點而言，它們的凸包就是由此三點所組成。但是接下來加入D時會發現，線段<K,D>；才會在凸包上，所以將線段<K,C>；排除，C點不可能是凸包。
⒊ 即當加入一點時，必須考慮到前面的線段是否在凸包上。從基點開始，凸包上每條相臨的線段的旋轉方向應該一致，並與掃描的方向相反。如果發現新加的點使得新線段與上線段的旋轉方向發生變化，則可判定上一點必然不在凸包上。實現時可用向量叉積進行判斷，設新加入的點為
，上一點為
，再上一點為
。順時針掃描時，如果向量
與
的叉積為正（逆時針掃描判斷是否為負），則將上一點刪除。刪除過程需要回溯，將之前所有叉積符號相反的點都刪除，然後將新點加入凸包。
在圖2中，加入K點時，由於線段<H,C>要旋轉到<H,K>的角度，為順時針旋轉，所以C點不在凸包上，應該刪除，保留K點。接着加入D點，由於線段<K,D>要旋轉到<H,K>的角度，為逆時針旋轉，故D點保留。按照上述步驟進行掃描，直到點集中所有的點都遍歷完成，即得到凸包。
這個算法可以直接在原數據上進行運算，因此空間複雜度為O⑴。但如果將凸包的結果存儲到另一數組中，則可能在代碼級別進行優化。由於在掃描凸包前要進行排序，因此時間複雜度至少為快速排序的O(nlgn)。後面的掃描過程複雜度為O(n)，因此整個算法的複雜度為O(nlgn)。
Jarvis步進法
對於一個有三個或以上點的點集Q，過程如下：
圖3 示例圖
計算點集最右邊的點為凸包的頂點的起點，如圖3的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];//保存凸包上的點，相當於所説的棧S
POINTa[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.1415827
typedefstructCseg
{
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; begin if abs(x-midx)<=zero then begin if abs(midy-y)<=zero then exit(0); if midy>y then exit(1) else exit(2); end else 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; var ax,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; var i,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; begin init; main; end.
