关闭

5.30 convexhull.cpp 学习

98人阅读 评论(0) 收藏 举报
分类:


转载:http://baike.baidu.com/link?url=LuzOfDLqjFSUfuT-hhug0NGKMKGdVTgNqzrGW8HMasR_nXbNDLmfqyYIGLmR7uP2U6jpeAV1P1d_6HdcGQSaqq

凸包

编辑 锁定
凸包(Convex Hull)是一个计算几何(图形学)中的概念。
在一个实数向量空间V中,对于给定集合X,所有包含X的凸集交集S被称为X的凸包
X的凸包可以用X内所有点(X1,...Xn)的线性组合来构造.
在二维欧几里得空间中,凸包可想象为一条刚好包著所有点的橡皮圈。
用不严谨的话来讲,给定二维平面上的点集,凸包就是将最外层的点连接起来构成的凸多边型,它能包含点集中所有的点。
中文名
凸包
外文名
Convex Hull
定    义
包含X的凸集交集
内    容
一个计算几何(图形学)中的概念
条    件
二维欧几里得空间

凸包定义

编辑
⒈对于一个集合D,D中任意有限个点的线性组合的全体称为D的凸包。
⒉对于一个集合D,所有包含D的凸集之交称为D的凸包。
可以证明,上述两种定义是等价的
概念
1  点集Q的凸包(convex hull)是指一个最小凸多边形,满足Q中的点或者在多边形边上或者在其内。右图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。
2  一组平面上的点,求一个包含所有点的最小的凸多边形,这就是凸包问题了。这可以形象地想成这样:在地上放置一些不可移动的木桩,用一根绳子把他们尽量紧地圈起来,并且为凸边形,这就是凸包了。

凸包平面求法

编辑

凸包常见求法

凸包最常用的凸包算法是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点不可能是凸包。
⒊ 即当加入一点时,必须考虑到前面的线段是否会出现在凸包上。从基点开始,凸包上每条相临的线段的旋转方向应该一致,并与扫描的方向相反。如果发现新加的点使得新线段与上线段的旋转方向发生变化,则可判定上一点必然不在凸包上。实现时可用向量叉积进行判断,设新加入的点为pn + 1,上一点为pn,再上一点为pn - 1。顺时针扫描时,如果向量<pn - 1,pn>;与<pn,pn + 1>;的叉积为正(逆时针扫描判断是否为负),则将上一点删除。删除过程需要回溯,将之前所有叉积符号相反的点都删除,然后将新点加入凸包。
在上图中,加入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)。[1] 

凸包代码实例

编辑

凸包代码一(C)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#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;
}
1
 

凸包代码二(C)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#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)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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.
0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:7606次
    • 积分:183
    • 等级:
    • 排名:千里之外
    • 原创:11篇
    • 转载:1篇
    • 译文:0篇
    • 评论:2条
    文章分类
    最新评论