POJ 3608
1.计算P上y坐标值最小的顶点(称为 yminP )和Q上y坐标值最大的顶点(称为 ymaxQ)。
2.为多边形在 yminP 和 ymaxQ 处构造两条切线 LP 和 LQ 使得他们对应的多边形位于他们的右侧。
此时 LP 和 LQ 拥有不同的方向, 并且 yminP 和 ymaxQ 成为了多边形间的一个对踵点对。
3.计算距离(yminP,ymaxQ) 并且将其维护为当前最小值。
4.顺时针同时旋转平行线直到其中一个与其所在的多边形的边重合。
5.如果只有一条线与边重合, 那么只需要计算“顶点-边”对踵点对和“顶点-顶点”对踵点对距离。 都将他们与当前最小值
比较, 如果小于当前最小值则进行替换更新。如果两条切线都与边重合,那么情况就更加复杂了。如果边“交叠”,也就是
可以构造一条与两条边都相交的公垂线(但不是在顶点处相交), 那么就计算“边-边”距离。 否则计算三个新的“顶点-顶
点”对踵点对距离。 所有的这些距离都与当前最小值进行比较, 若小于当前最小值则更新替换。
6.重复执行步骤4和步骤5, 直到新的点对为(yminP,ymaxQ)。
7.输出最小距离
这是求两凸包最短距离的经典算法。但是,不知为什么,我的代码死活过不了。。。T_T
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <algorithm> 5 #include <cmath> 6 using namespace std; 7 const double eps=0.00000001; 8 9 struct point { 10 double x,y; 11 }p[10050],q[10050]; 12 int n,m; 13 int ansp[10050],ansq[10050],cntp,cntq; 14 int st[10050],stop; 15 16 bool cmp(point A,point B){ 17 if(A.y<B.y) return true; 18 else if(A.y==B.y){ 19 if(A.x<B.x) return true; 20 } 21 return false; 22 } 23 24 double dist(point a , point b){ 25 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); 26 } 27 28 double multi(point a, point b, point c){ 29 point p1; p1.x=a.x-c.x; p1.y=a.y-c.y; 30 point p2; p2.x=b.x-c.x; p2.y=b.y-c.y; 31 return (p1.x*p2.y-p1.y*p2.x); 32 } 33 34 void forTU(point *pot, int &allVex, int *anp, int &cnt){ 35 // cout<<allVex<<"all"<<endl; 36 stop=cnt=0; 37 st[stop++]=0; st[stop++]=1; 38 for(int i=2;i<allVex;i++){ 39 while(stop>1&&multi(pot[i],pot[st[stop-1]],pot[st[stop-2]])<=0) stop--; 40 st[stop++]=i; 41 } 42 for(int i=0;i<stop;i++) 43 anp[cnt++]=st[i]; 44 stop=0; 45 st[stop++]=allVex-1; st[stop++]=allVex-2; 46 for(int i=allVex-3;i>=0;i--){ 47 while(stop>1&&multi(pot[i],pot[st[stop-1]],pot[st[stop-2]])<=0) stop--; 48 st[stop++]=i; 49 } 50 for(int i=1;i<stop;i++) 51 anp[cnt++]=st[i]; 52 // for(int i=0;i<cnt;i++) 53 // cout<<anp[i]<<endl; 54 // cout<<endl; 55 } 56 57 double dotcross(point a,point b, point c){ 58 point p1; p1.x=a.x-c.x; p1.y=a.y-c.y; 59 point p2; p2.x=b.x-c.x; p2.y=b.y-c.y; 60 return p1.x*p2.x+p1.y*p2.y; 61 } 62 63 double pline(point a,point b,point c){ 64 return (fabs(multi(a,b,c)))/(dist(a,b)); 65 } 66 67 double pseg(point a,point b,point c){ 68 if(dotcross(a,c,b)<-eps) return dist(b,c); 69 if(dotcross(b,c,a)<-eps) return dist(a,c); 70 return pline(a,b,c); 71 } 72 73 double paral(point a1,point a2, point b1,point b2 ){ 74 double ans1=min(pseg(a1,a2,b1),pseg(a1,a2,b2)); 75 double ans2=min(pseg(b1,b2,a1),pseg(b1,b2,a2)); 76 return min(ans1,ans2); 77 } 78 79 int sgn(double x) 80 { 81 if(fabs(x) < eps)return 0; 82 if(x < 0)return -1; 83 else return 1; 84 } 85 86 double Get_angle(point a1,point a2,point b1,point b2) 87 { 88 point p1; p1.x=a1.x-a2.x; p1.y=a1.y-a2.y; 89 point p2; p2.x=b1.x-b2.x; p2.y=b1.y-b2.y; 90 return p1.x*p2.y-p1.y*p2.x; 91 } 92 93 double slove(point *p, int *ap, int &cp, point *q, int *aq, int &cq){ 94 int sp=0,sq=0; double tmp; 95 for(int i=0;i<cq;i++) //max 96 if(q[aq[i]].y-eps>q[aq[sq]].y) sq=i; 97 double res=dist(p[ap[sp]],q[aq[sq]]); 98 for(int i=0;i<cp;i++){ 99 // while((tmp=fabs(multi(p[ap[i]],p[ap[i+1]],q[aq[sq]])/2)-fabs(multi(p[ap[i]],p[ap[i+1]],q[aq[(sq+1)%cq]])/2))>eps) 100 // while(tmp=multi(p[ap[i]],p[ap[i+1]],q[aq[(sq+1)%cq]])-multi(p[ap[i]],p[ap[i+1]],q[aq[sq]])>eps) 101 // sq=(sq+1)%cq; 102 // if(tmp<-eps){ 103 while(sgn(tmp = Get_angle(p[i],p[(i+1)%cp],q[sq],q[(sq+1)%cq])) < 0 ) 104 sq = (sq + 1)%cq; 105 if(sgn(tmp) == 0) 106 res=min(res,pseg(p[ap[i]],p[ap[i+1]],q[aq[sq]])); 107 // cout<<res<<endl; 108 // } 109 else{ 110 res=min(res,paral(p[ap[i]],p[ap[i+1]],q[aq[sq]],q[aq[(sq+1)%cq]])); 111 // cout<<res<<endl; 112 } 113 } 114 return res; 115 } 116 117 118 int main(){ 119 double ans=10000000; 120 while(scanf("%d%d",&n,&m)!=EOF){ 121 if(n==0&&m==0) break; 122 for(int i=0;i<n;i++) 123 scanf("%lf%lf",&p[i].x,&p[i].y); 124 for(int i=0;i<m;i++) 125 scanf("%lf%lf",&q[i].x,&q[i].y); 126 sort(p,p+n,cmp); 127 sort(q,q+m,cmp); 128 forTU(p,n,ansp,cntp); 129 forTU(q,m,ansq,cntq); 130 ans=1e99; 131 ans=min(ans,slove(p,ansp,cntp,q,ansq,cntq)); //min,max 132 ans=min(ans,slove(q,ansq,cntq,p,ansp,cntp)); 133 printf("%.5lf\n",ans); 134 } 135 return 0; 136 }
别人的代码:
#include <iostream> #include <string.h> #include <algorithm> #include <stdio.h> #include <math.h> using namespace std; const int N=50000; const double eps=1e-9; const double INF=1e99; struct Point { double x,y; }; Point P[N],Q[N]; double cross(Point A,Point B,Point C) { return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x); } double dist(Point A,Point B) { return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)); } double multi(Point A,Point B,Point C) { return (B.x-A.x)*(C.x-A.x)+(B.y-A.y)*(C.y-A.y); } //顺时针排序 void anticlockwise(Point p[],int n) { for(int i=0;i<n-2;i++) { double tmp=cross(p[i],p[i+1],p[i+2]); if(tmp>eps) return; else if(tmp<-eps) { reverse(p,p+n); return; } } } //计算C点到直线AB的最短距离 double Getdist(Point A,Point B,Point C) { if(dist(A,B)<eps) return dist(B,C); if(multi(A,B,C)<-eps) return dist(A,C); if(multi(B,A,C)<-eps) return dist(B,C); return fabs(cross(A,B,C)/dist(A,B)); } //求一条直线的两端点到另外一条直线的距离,反过来一样,共4种情况 double MinDist(Point A,Point B,Point C,Point D) { return min(min(Getdist(A,B,C),Getdist(A,B,D)),min(Getdist(C,D,A),Getdist(C,D,B))); } double Solve(Point P[],Point Q[],int n,int m) { int yminP=0,ymaxQ=0; for(int i=0;i<n;i++) if(P[i].y<P[yminP].y) yminP=i; for(int i=0;i<m;i++) if(Q[i].y>Q[ymaxQ].y) ymaxQ=i; P[n]=P[0]; Q[m]=Q[0]; double tmp,ans=INF; for(int i=0;i<n;i++) { while(tmp=cross(P[yminP+1],Q[ymaxQ+1],P[yminP])-cross(P[yminP+1],Q[ymaxQ],P[yminP])>eps) ymaxQ=(ymaxQ+1)%m; if(tmp<-eps) ans=min(ans,Getdist(P[yminP],P[yminP+1],Q[ymaxQ])); else ans=min(ans,MinDist(P[yminP],P[yminP+1],Q[ymaxQ],Q[ymaxQ+1])); yminP=(yminP+1)%n; } return ans; } int main() { int n,m; while(cin>>n>>m) { if(n==0&&m==0) break; for(int i=0;i<n;i++) cin>>P[i].x>>P[i].y; for(int i=0;i<m;i++) cin>>Q[i].x>>Q[i].y; anticlockwise(P,n); anticlockwise(Q,m); printf("%.5lf\n",min(Solve(P,Q,n,m),Solve(Q,P,m,n))); } return 0; }