本文共 6003 字,大约阅读时间需要 20 分钟。
此文章是一些自己未接触过的计算几何知识,不是太明白。。。。。
&& 本题代码来自 大牛 :
/** 最小圆覆盖 (随机增量法) 求给你n个点,让你用一个最小的圆覆盖所有的点,求这个圆的坐标&&半径**/#include#include #include #include #include using namespace std;const double eps = 1e-12;const int len = 100010;struct Point{ double x,y;} p[len];struct Line{ Point a,b;};int dbcmp(double n){ return n < -eps ? -1 : n > eps;}double dis(Point a, Point b){ return sqrt((a.x-b.x) * ( a.x-b.x) + ( a.y-b.y) * ( a.y-b.y));}//求两直线的交点Point intersection(Line u,Line v){ Point ret=u.a; double t=((u.a.x-v.a.x)*(v.b.y-v.a.y)-(u.a.y-v.a.y)*(v.b.x-v.a.x)) /((u.a.x-u.b.x)*(v.b.y-v.a.y)-(u.a.y-u.b.y)*(v.b.x-v.a.x)); ret.x+=(u.b.x-u.a.x)*t; ret.y+=(u.b.y-u.a.y)*t; return ret;}//三角形外接圆圆心(外心)Point center(Point a,Point b,Point c){ Line u,v; u.a.x=(a.x+b.x)/2; u.a.y=(a.y+b.y)/2; u.b.x=u.a.x+(u.a.y-a.y); u.b.y=u.a.y-(u.a.x-a.x); v.a.x=(a.x+c.x)/2; v.a.y=(a.y+c.y)/2; v.b.x=v.a.x+(v.a.y-a.y); v.b.y=v.a.y-(v.a.x-a.x); return intersection(u,v);}void min_cir(Point * p, int n, Point & cir, double &r){ random_shuffle(p, p + n); cir = p[0]; r = 0; for(int i = 1; i < n; ++i) { if(dbcmp(dis(p[i],cir) -r) <= 0)continue; cir = p[i]; r = 0; for(int j =0; j < i ; ++j) { if(dbcmp(dis(p[j], cir) -r) <= 0 )continue; cir.x = (p[i].x + p[j].x) /2; cir.y = (p[i].y + p[j].y) /2; r = dis( cir, p[j]); for(int k =0; k < j; ++k) { if(dbcmp( dis(p[k], cir) -r) <= 0) continue; cir = center(p[i], p[j], p[k]); r = dis( cir, p[k]); } } }}int main(){ int n; while( ~scanf("%d", &n), n) { for (int i = 0; i < n; ++i) { scanf("%lf%lf", &p[i].x, &p[i].y); } Point cir; double r; min_cir(p, n, cir, r); printf("%.2lf %.2lf %.2lf\n", cir.x, cir.y, r); }}
题意:求两凸包之间的距离;摘自 kuangbin 大牛:
#include#include #include #include #include using namespace std;const double eps = 1e-8;int sgn(double x){ if(fabs(x) < eps)return 0; if(x < 0)return -1; else return 1;}struct Point{ double x,y; Point(double _x = 0.0,double _y = 0.0){ x = _x;y = _y; } Point operator -(const Point &b)const{ return Point(x - b.x, y - b.y); } double operator ^(const Point &b)const{ return x*b.y - y*b.x; } double operator *(const Point &b)const{ return x*b.x + y*b.y; } void input(){ scanf("%lf%lf",&x,&y); }};struct Line{ Point s,e; Line(){} Line(Point _s,Point _e){ s = _s; e = _e; }};double dist(Point a,Point b){ return sqrt((a-b)*(a-b));}Point NearestPointToLineSeg(Point P,Line L){ Point result; double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s)); if(t >= 0 && t <= 1){ result.x = L.s.x + (L.e.x - L.s.x)*t; result.y = L.s.y + (L.e.y - L.s.y)*t; } else{ if(dist(P,L.s) < dist(P,L.e)) result = L.s; else result = L.e; } return result;}/* * 求凸包,Graham算法 * 点的编号0~n-1 * 返回凸包结果Stack[0~top-1]为凸包的编号 */const int MAXN = 10010;Point list[MAXN];int Stack[MAXN],top;//相对于list[0]的极角排序bool _cmp(Point p1,Point p2){ double tmp = (p1-list[0])^(p2-list[0]); if(sgn(tmp) > 0)return true; else if(sgn(tmp) == 0 && sgn(dist(p1,list[0]) - dist(p2,list[0])) <= 0) return true; else return false;}void Graham(int n){ Point p0; int k = 0; p0 = list[0]; //找最下边的一个点 for(int i = 1;i < n;i++){ if( (p0.y > list[i].y) || (p0.y == list[i].y && p0.x > list[i].x) ){ p0 = list[i]; k = i; } } swap(list[k],list[0]); sort(list+1,list+n,_cmp); if(n == 1){ top = 1; Stack[0] = 0; return; } if(n == 2){ top = 2; Stack[0] = 0; Stack[1] = 1; return ; } Stack[0] = 0; Stack[1] = 1; top = 2; for(int i = 2;i < n;i++){ while(top > 1 && sgn((list[Stack[top-1]]-list[Stack[top-2]])^(list[i]-list[Stack[top-2]])) <= 0) top--; Stack[top++] = i; }}//点p0到线段p1p2的距离double pointtoseg(Point p0,Point p1,Point p2){ return dist(p0,NearestPointToLineSeg(p0,Line(p1,p2)));}//平行线段p0p1和p2p3的距离double dispallseg(Point p0,Point p1,Point p2,Point p3){ double ans1 = min(pointtoseg(p0,p2,p3),pointtoseg(p1,p2,p3)); double ans2 = min(pointtoseg(p2,p0,p1),pointtoseg(p3,p0,p1)); return min(ans1,ans2);}//得到向量a1a2和b1b2的位置关系double Get_angle(Point a1,Point a2,Point b1,Point b2){ Point t = b1 - ( b2 - a1 ); return (a2-a1)^(t-a1);}//旋转卡壳,求两个凸包的最小距离double rotating_calipers(Point p[],int np,Point q[],int nq){ int sp = 0, sq = 0; for(int i = 0;i < np;i++) if(sgn(p[i].y - p[sp].y) < 0) sp = i; for(int i = 0;i < nq;i++) if(sgn(q[i].y - q[sq].y) > 0) sq = i; double tmp; double ans = 1e99; for(int i = 0;i < np;i++){ while(sgn(tmp = Get_angle(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq])) < 0 ) sq = (sq + 1)%nq; if(sgn(tmp) == 0) ans = min(ans,dispallseg(p[sp],p[(sp+1)%np],q[sq],q[(sq+1)%nq])); else ans = min(ans,pointtoseg(q[sq],p[sp],p[(sp+1)%np])); sp = (sp+1)%np; } return ans;}double solve(Point p[],int n,Point q[],int m){ return min(rotating_calipers(p,n,q,m),rotating_calipers(q,m,p,n));}Point p[MAXN],q[MAXN];int main(){ int n,m; while(scanf("%d%d",&n,&m)==2) { if(n == 0 && m == 0)break; for(int i = 0;i < n;i++) list[i].input(); Graham(n); n = top; for(int i = 0;i < n;i++) p[i] = list[Stack[i]]; for(int i = 0;i < m;i++) list[i].input(); Graham(m); m = top; for(int i = 0;i < m;i++) q[i] = list[Stack[i]]; printf("%.5lf\n",solve(p,n,q,m)); } return 0;}
转载地址:http://qcsgi.baihongyu.com/