凸包就是把给定点包围在内部,面积最小的凸多边形
Graham扫描法
1.把点按照x坐标排序
2.将p1,p2放到凸包中
3.从p3开始,当新的点在凸包的前进方向的左边时加入该点,否则弹出最后加入的点,直到新点在左边
4.将上述流程反过来做一遍求上凸壳
struct vector { double x,y; vector (double X=0,double Y=0) { x=X,y=Y; } } typedef vector point; double cross(vector a,vector b) { return a.x*b.y-a.y*b.x; } void convexhull(point *p,int n,point *ch) { sort(p+1,p+n+1); if (n==1){ ch[++m]=n; return ; } m=1; for (int i=1;i<=n;i++){ while (m>2&&dcmp(cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<=0) m--; ch[m++]=p[i]; } int k=m; for (int i=n-1;i>=1;i--){ while (m>k&&dcmp(cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<=0) m--; ch[m++]=p[i]; } m--; }
旋转卡壳
求凸包上最远点对的距离
被凸包上一对平行直线卡住的点对叫做对踵点,答案一定在对踵点上。
按照逆时针顺序枚举凸包上一条边,则凸包上到该边所在直线最远的点也单调逆时针旋转,相当于用一条直线卡对面一个顶点
1 double rotating_calipers(point* ch,int n)//旋转卡壳,被凸包上被一对平行直线卡住的点叫对踵点,最远点对一定在凸包的一对对踵点上 2 { 3 if (n==1) return 0; 4 if (n==2) return dot(ch[0]-ch[1],ch[0]-ch[1]);//如果只有两个点,那么就是两点的直线距离 5 int now=1,i; 6 double ans=0; 7 ch[n]=ch[0]; 8 for (int i=0;i<n;i++) 9 { 10 while(dcmp(distl(ch[now],ch[i],ch[i+1])-distl(ch[now+1],ch[i],ch[i+1]))<=0)//最远点随着平行线的旋转是单调的,所以点不会来回移动 11 now=(now+1)%n; 12 ans=max(ans,dot(ch[now]-ch[i],ch[now]-ch[i])); 13 ans=max(ans,dot(ch[now]-ch[i+1],ch[now]-ch[i+1]));//找到与当前直线平行的最远点,用来更新答案 14 } 15 return ans; 16 }
求两凸包上的最近距离
1 double rotating_calipers(point* p,point* q,int n,int m)//利用旋转卡壳求两个凸包上最近点的距离,一个凸包上的平行线逆时针旋转,另一个凸包上的最远点也单调逆时针旋转,所以这个算法要正反进行两遍 2 { 3 int x=0; int y=0,i; 4 double ans=1e10,t; 5 for (int i=0;i<n;i++) 6 x=(p[i].y<p[x].y)?i:x; 7 for (int i=0;i<m;i++) 8 y=(p[i].y>p[y].y)?i:y; 9 p[n]=p[0]; q[m]=q[0]; 10 for (int i=0;i<n;i++) 11 { 12 while((t=dcmp(cross(p[x]-p[x+1],q[y+1]-q[y])))<0)//判断凸壳上下一条边的旋转方向 13 y=(y+1)%m; 14 if (t==0)//平行的时候四个点之间更新答案 15 { 16 ans=min(ans,dists(p[x],q[y+1],q[y])); 17 ans=min(ans,dists(p[x+1],q[y+1],q[y])); 18 ans=min(ans,dists(q[y],p[x],p[x+1])); 19 ans=min(ans,dists(q[y+1],p[x],p[x+1])); 20 } 21 else 22 ans=min(ans,dists(q[y],p[x],p[x+1]));//否则只用最远点更新答案 23 x=(x+1)%n; 24 } 25 return ans; 26 }
最小圆覆盖
用最小的圆覆盖平面中的所有点。
1 point center(point a,point b,point c){ 2 point p=(a+b)/2; point q=(a+c)/2; 3 vector v=rotate(b-a,pi/2); vector u=rotate(c-a,pi/2); 4 if (dcmp(cross(v,u))==0) { 5 if(dcmp(len(a-b)+len(b-c)-len(a-c))==0) 6 return (a+c)/2; 7 if(dcmp(len(a-c)+len(b-c)-len(a-b))==0) 8 return (a+b)/2; 9 if(dcmp(len(a-b)+len(a-c)-len(b-c))==0) 10 return (b+c)/2; 11 } 12 return glt(p,v,q,u); 13 } 14 double mincc(point *p,int n,point &c) 15 { 16 random_shuffle(p,p+n); 17 c=p[0]; 18 double r=0; 19 int i,j,k; 20 for (i=1;i<n;i++) 21 if (dcmp(len(c-p[i])-r)>0){ 22 c=p[i],r=0; 23 for (j=0;j<i;j++) 24 if (dcmp(len(c-p[j])-r)>0){ 25 c=(p[i]+p[j])/2; 26 r=len(c-p[i]); 27 for (k=0;k<j;k++) 28 if (dcmp(len(c-p[k])-r)>0){ 29 c=center(p[i],p[j],p[k]); 30 r=len(c-p[i]); 31 } 32 } 33 } 34 return r; 35 }