寻找凸包的graham 扫描法

1,点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。 
2,凸包最常用的凸包算法是Graham扫描法和Jarvis步进法。 
3,Graham扫描法： 
首先，找到所有点中最左边的(y坐标最小的)，如果y坐标相同，找x坐标最小的. 
以这个点为基准求所有点的极角（atan2(y-y0,x-x0)），并按照极角对这些点排序，前述基准点在最前面，设这些点为P[0]..P[n-1. 
注:这样预处理后,保证p[0],p[1]和p[n-1]都是凸包上的点. 
建立一个栈,初始时P[0]、P[1]、P[2]进栈,对于 P[3..n-1]的每个点，若栈顶的两个点与它不构成"向左转"的关系,则将栈顶的点出栈,直至没有点需要出栈以后将当前点进栈； 
所有点处理完之后栈中保存的点就是凸包了。 
图示: 

 

4,实现代码 

Cpp代码  

1. #include <iostream>  
2. #include <cmath>  
3. using namespace std;  
4.   
5. /* 
6. PointSet[]：输入的点集 
7. ch[]：输出的凸包上的点集，按照逆时针方向排列 
8. n：PointSet中的点的数目 
9. len：输出的凸包上的点的个数 
10. */  
11.   
12. struct Point  
13. {  
14.     float x,y;  
15. };  
16.   
17. //小于0,说明向量p0p1的极角大于p0p2的极角  
18. float multiply(Point p1,Point p2,Point p0)  
19. {  
20.     return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));  
21. }  
22.   
23. float dis(Point p1,Point p2)  
24. {  
25.     return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));  
26. }  
27.   
28.   
29. void Graham_scan(Point PointSet[],Point ch[],int n,int &len)  
30. {  
31.     int i,j,k=0,top=2;  
32.     Point tmp;  
33.   
34.     //找到最下且偏左的那个点  
35.     for(i=1;i<n;i++)  
36.         if ((PointSet[i].y<PointSet[k].y)||((PointSet[i].y==PointSet[k].y)&&(PointSet[i].x<PointSet[k].x)))  
37.             k=i;  
38.     //将这个点指定为PointSet[0]  
39.     tmp=PointSet[0];  
40.     PointSet[0]=PointSet[k];  
41.     PointSet[k]=tmp;  
42.   
43.     //按极角从小到大,距离偏短进行排序  
44.     for (i=1;i<n-1;i++)  
45.     {  
46.         k=i;  
47.         for (j=i+1;j<n;j++)  
48.             if( (multiply(PointSet[j],PointSet[k],PointSet[0])>0)  
49.                 ||((multiply(PointSet[j],PointSet[k],PointSet[0])==0)  
50.                     &&(dis(PointSet[0],PointSet[j])<dis(PointSet[0],PointSet[k]))) )  
51.                 k=j;//k保存极角最小的那个点,或者相同距离原点最近  
52.         tmp=PointSet[i];  
53.         PointSet[i]=PointSet[k];  
54.         PointSet[k]=tmp;  
55.     }  
56.     //第三个点先入栈  
57.     ch[0]=PointSet[0];  
58.     ch[1]=PointSet[1];  
59.     ch[2]=PointSet[2];  
60.     //判断与其余所有点的关系  
61.     for (i=3;i<n;i++)  
62.     {  
63.         //不满足向左转的关系,栈顶元素出栈  
64.         while(multiply(PointSet[i],ch[top],ch[top-1])>=0) top--;  
65.         //当前点与栈内所有点满足向左关系,因此入栈.  
66.         ch[++top]=PointSet[i];  
67.     }  
68.     len=top+1;  
69. }  
70.   
71. const int maxN=1000;  
72. Point PointSet[maxN];  
73. Point ch[maxN];  
74. int n;  
75. int len;  
76.   
77. int main()  
78. {  
79.     int n=5;  
80.     float x[]={0,3,4,2,1};  
81.     float y[]={0,0,0,3,1};  
82.     for(int i=0;i<n;i++)  
83.     {  
84.         PointSet[i].x=x[i];  
85.         PointSet[i].y=y[i];  
86.     }  
87.     Graham_scan(PointSet,ch,n,len);  
88.     for(int i=0;i<len;i++)  
89.         cout<<ch[i].x<<" "<<ch[i].y<<endl;  
90.     return 0;  
91. }  

5,应用: http://acm.pku.edu.cn/JudgeOnline/problem?id=1113  
实现代码: 

Cpp代码  

1. #include <iostream>  
2. #include <cmath>  
3. using namespace std;  
4.   
5. /* 
6. PointSet[]：输入的点集 
7. ch[]：输出的凸包上的点集，按照逆时针方向排列 
8. n：PointSet中的点的数目 
9. len：输出的凸包上的点的个数 
10. */  
11.   
12. struct Point  
13. {  
14.     float x,y;  
15. };  
16.   
17. //小于0,说明向量p0p1的极角大于p0p2的极角  
18. float multiply(Point p1,Point p2,Point p0)  
19. {  
20.     return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));  
21. }  
22.   
23. float dis(Point p1,Point p2)  
24. {  
25.     return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));  
26. }  
27.   
28.   
29. void Graham_scan(Point PointSet[],Point ch[],int n,int &len)  
30. {  
31.     int i,j,k=0,top=2;  
32.     Point tmp;  
33.   
34.     //找到最下且偏左的那个点  
35.     for(i=1;i<n;i++)  
36.         if ((PointSet[i].y<PointSet[k].y)||((PointSet[i].y==PointSet[k].y)&&(PointSet[i].x<PointSet[k].x)))  
37.             k=i;  
38.     //将这个点指定为PointSet[0]  
39.     tmp=PointSet[0];  
40.     PointSet[0]=PointSet[k];  
41.     PointSet[k]=tmp;  
42.   
43.     //按极角从小到大,距离偏短进行排序  
44.     for (i=1;i<n-1;i++)  
45.     {  
46.         k=i;  
47.         for (j=i+1;j<n;j++)  
48.             if( (multiply(PointSet[j],PointSet[k],PointSet[0])>0)  
49.                 ||((multiply(PointSet[j],PointSet[k],PointSet[0])==0)  
50.                     &&(dis(PointSet[0],PointSet[j])<dis(PointSet[0],PointSet[k]))) )  
51.                 k=j;//k保存极角最小的那个点,或者相同距离原点最近  
52.         tmp=PointSet[i];  
53.         PointSet[i]=PointSet[k];  
54.         PointSet[k]=tmp;  
55.     }  
56.     //第三个点先入栈  
57.     ch[0]=PointSet[0];  
58.     ch[1]=PointSet[1];  
59.     ch[2]=PointSet[2];  
60.     //判断与其余所有点的关系  
61.     for (i=3;i<n;i++)  
62.     {  
63.         //不满足向左转的关系,栈顶元素出栈  
64.         while(multiply(PointSet[i],ch[top],ch[top-1])>=0) top--;  
65.         //当前点与栈内所有点满足向左关系,因此入栈.  
66.         ch[++top]=PointSet[i];  
67.     }  
68.     len=top+1;  
69. }  
70.   
71. const int maxN=1010;  
72. Point PointSet[maxN];  
73. Point ch[maxN];  
74. int n;  
75. int len;  
76.   
77. int main()  
78. {  
79.     double ans=0;  
80.     int d;  
81.     cin>>n>>d;  
82.     for(int i=0;i<n;i++)  
83.         cin>>PointSet[i].x>>PointSet[i].y;//input the data;  
84.     Graham_scan(PointSet,ch,n,len);  
85.     for(int i=0;i<len;i++)  
86.         ans+=dis(ch[i],ch[(i+1)%len]);  
87.     ans+=2*d*acos(-1.0); //等价于圆形周长  
88.     cout<<(int)(ans+0.5)<<endl; //四舍五入  
89.     return 0;  
90. }  

附：矢量叉积

   设有点p0(x0,y0),p1(x1,y1),p2(x2,y2).（p0p1）,（p0p2）是共p0的两条向量，叉积d = （p0p1）x(p0p2) = (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0)

   叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系：
若 d > 0 , 则(p0p1)在(p0p2)的顺时针方向。 
若 d < 0 , 则(p0p1)在(p0p2)的逆时针方向。(图示方向)
若 d = 0 , 则(p0p1)与(p0p2)共线，但可能同向也可能反向。

                                        

FROM: http://kmplayer.iteye.com/blog/604405