凸包

概念1.1 点集Q的凸包(convex hull)是指一个最小凸多边形,满足Q中的点或者在多边形边上或者在其内。右图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。
1.2 一组平面上的点,求一个包含所有点的最小的凸多边形,这就是凸包问题了。这可以形象地想成这样:在地上放置一些不可移动的木桩,用一根绳子把他们尽量紧地圈起来,这就是凸包了。
平面凸包的求法2.1凸包最常用的凸包算法是Graham扫描法和Jarvis步进法。
对于一个有三个或以上点的点集Q,过程如下:
计算点集最右边的点为凸包的顶点的起点,如上图的P3点。
Do
For i = 0 To 总顶点数
计算有向向量P3->Pi
If 其余顶点全部在有向向量P3->Pi的左侧或右侧,则Pi点为凸包的下一顶点
Pi点加入凸包列表
GoTo 1
End If
Next
Exit Do
1:
Loop
此过程执行后,点按极角自动顺时针或逆时针排序,只需要按任意两点的次序就可以了。而左侧或右侧的判断可以用前述的矢量点积性质实现。
2.2求凸包有很多方法,不过最适合OI的估计还是Graham's Scan这个方法了。它的大致方法是这样的: 首先,找到所有点中最左边的(y坐标最小的),如果y坐标相同,找x坐标最小的;以这个点为基准求所有点的极角(atan2(y-y0,x-x0)),并按照极角对这些点排序,前述基准点在最前面,设这些点为P[0]..P[n-1];建立一个栈,初始时P[0]、P[1]、P[2]进栈,对于P[3..n-1]的每个点,若栈顶的两个点与它不构成“向左转”的关系,则将栈顶的点出栈,直至没有点需要出栈以后将当前点进栈;所有点处理完之后栈中保存的点就是凸包了。
如何判断A、B、C构成的关系不是向左转呢?如果b-a与c-a的叉乘小于0就不是。a与b的叉乘就是a.x*b.y-a.y*b.x。
上面的这个Graham的实现比我原来按照USACO里的课文写得简单多了,主要是它通过简单的预处理保证了P[0]、P[1]以及P[n-1]肯定是凸包里的点,这样就可以避免在凸包“绕回来”的时候繁杂的处理
代码===============C++代码===================
#include <iostream> // 求点集合的凸包的gram算法。n是顶点个数,x,y是顶点
坐标。
#include <fstream> // order 是按照顶点和左下脚的角度的排序后数组。
#include <deque> // tu即是逆时针的凸包上的顶点。
#include <math.h> //
using namespace std; //使用条件: 1。点可以任意给,可重复。
// 2。三个以及以上的点。
ifstream fin("gram.in"); // 3。已经考虑了边上有点的情况。
#define NN 1000
#define pi 3.1415827
typedef struct Cseg{
double x,y,tg;
}Cseg;
int n;
double x[NN],y[NN];
deque <Cseg> order;
deque <int> tu;
Cseg seg1;
deque <Cseg> ::iterator p1;
deque <int> ::iterator p,q;
void in();
void gram();
void makeorder(int s);
double dist(double x1,double yy1,double x2,double yy2);
double cross(double x1,double yy1,double x2,double yy2);
void out();
int main()
{
in();
gram();
out();
return 0;
}
void out()
{
int i;
for (i=0;i<tu.size();i++){
cout<<order[tu].x<<" "<<order[tu].y<<endl;
}
cout<<tu.size()<<" Edges Polydgon"<<endl;
return;
}
void in()
{
int i;
fin>>n;
for (i=0;i<n;i++)
fin>>x>>y;
return;
}
void gram()
{
int i,mm;
mm=0;
for (i=1;i<n;i++)
if (y[mm]>y+1e-9) mm=i;
else if (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 (i=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;
//tu.push_back(i);
}
}//for
return;
}
void makeorder(int s)
{
int i;
double tg;
order.clear();
for (i=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;
}
else
if (tg<p1->tg){
order.insert(p1,seg1);
break;
}
p1++;
}//while
if (p1==order.end()) order.insert(p1,seg1);
}//for
seg1.x=x[s];seg1.y=y[s];
order.insert(order.begin(),seg1);
//for (i=0;i<order.size();i++)
// printf("i=%d %lf %lf
%lf
",i,order.x,order.y,order.tg*180/pi);
return;
}
double cross(double x1,double yy1,double x2,double yy2)
{
return (x1*yy2-x2*yy1);
}
double dist(double x1,double yy1,double x2,double yy2)
{
return pow((x1-x2)*(x1-x2)+(yy1-yy2)*(yy1-yy2),0.5);
}
======================================
P标程{pku 1113 }
{$Q-,S-,R-}
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,temp);
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)<-zero 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:0);
end;
begin
init;
main;
end.