poj3608 Bridge Across Islands
地址:http://poj.org/problem?id=3608
题目:
Time Limit: 1000MS | Memory Limit: 65536K | |||
Total Submissions: 11259 | Accepted: 3307 | Special Judge |
Description
Thousands of thousands years ago there was a small kingdom located in the middle of the Pacific Ocean. The territory of the kingdom consists two separated islands. Due to the impact of the ocean current, the shapes of both the islands became convex polygons. The king of the kingdom wanted to establish a bridge to connect the two islands. To minimize the cost, the king asked you, the bishop, to find the minimal distance between the boundaries of the two islands.
![](http://poj.org/images/3608_1.gif)
Input
The input consists of several test cases.
Each test case begins with two integers N, M. (3 ≤ N, M ≤ 10000)
Each of the next N lines contains a pair of coordinates, which describes the position of a vertex in one convex polygon.
Each of the next M lines contains a pair of coordinates, which describes the position of a vertex in the other convex polygon.
A line with N = M = 0 indicates the end of input.
The coordinates are within the range [-10000, 10000].
Output
For each test case output the minimal distance. An error within 0.001 is acceptable.
Sample Input
4 4 0.00000 0.00000 0.00000 1.00000 1.00000 1.00000 1.00000 0.00000 2.00000 0.00000 2.00000 1.00000 3.00000 1.00000 3.00000 0.00000 0 0
Sample Output
1.00000
Source
思路:
两凸包间的最近点对,套模板。
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <algorithm> 5 6 7 using namespace std; 8 const double PI = acos(-1.0); 9 const double eps = 1e-10; 10 11 /****************常用函数***************/ 12 //判断ta与tb的大小关系 13 int sgn( double ta, double tb) 14 { 15 if(fabs(ta-tb)<eps)return 0; 16 if(ta<tb) return -1; 17 return 1; 18 } 19 20 //点 21 class Point 22 { 23 public: 24 25 double x, y; 26 27 Point(){} 28 Point( double tx, double ty){ x = tx, y = ty;} 29 30 bool operator < (const Point &_se) const 31 { 32 return x<_se.x || (x==_se.x && y<_se.y); 33 } 34 friend Point operator + (const Point &_st,const Point &_se) 35 { 36 return Point(_st.x + _se.x, _st.y + _se.y); 37 } 38 friend Point operator - (const Point &_st,const Point &_se) 39 { 40 return Point(_st.x - _se.x, _st.y - _se.y); 41 } 42 //点位置相同(double类型) 43 bool operator == (const Point &_off)const 44 { 45 return sgn(x, _off.x) == 0 && sgn(y, _off.y) == 0; 46 } 47 48 }; 49 50 /****************常用函数***************/ 51 //点乘 52 double dot(const Point &po,const Point &ps,const Point &pe) 53 { 54 return (ps.x - po.x) * (pe.x - po.x) + (ps.y - po.y) * (pe.y - po.y); 55 } 56 //叉乘 57 double xmult(const Point &po,const Point &ps,const Point &pe) 58 { 59 return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.y); 60 } 61 //两点间距离的平方 62 double getdis2(const Point &st,const Point &se) 63 { 64 return (st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y); 65 } 66 //两点间距离 67 double getdis(const Point &st,const Point &se) 68 { 69 return sqrt((st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y)); 70 } 71 72 //两点表示的向量 73 class Line 74 { 75 public: 76 77 Point s, e;//两点表示,起点[s],终点[e] 78 double a, b, c;//一般式,ax+by+c=0 79 double angle;//向量的角度,[-pi,pi] 80 81 Line(){} 82 Line( Point ts, Point te):s(ts),e(te){}//get_angle();} 83 Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){} 84 85 //排序用 86 bool operator < (const Line &ta)const 87 { 88 return angle<ta.angle; 89 } 90 //向量与向量的叉乘 91 friend double operator / ( const Line &_st, const Line &_se) 92 { 93 return (_st.e.x - _st.s.x) * (_se.e.y - _se.s.y) - (_st.e.y - _st.s.y) * (_se.e.x - _se.s.x); 94 } 95 //向量间的点乘 96 friend double operator *( const Line &_st, const Line &_se) 97 { 98 return (_st.e.x - _st.s.x) * (_se.e.x - _se.s.x) - (_st.e.y - _st.s.y) * (_se.e.y - _se.s.y); 99 } 100 //从两点表示转换为一般表示 101 //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2 102 bool pton() 103 { 104 a = e.y - s.y; 105 b = s.x - e.x; 106 c = e.x * s.y - e.y * s.x; 107 return true; 108 } 109 //半平面交用 110 //点在向量左边(右边的小于号改成大于号即可,在对应直线上则加上=号) 111 friend bool operator < (const Point &_Off, const Line &_Ori) 112 { 113 return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x) 114 < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x); 115 } 116 //求直线或向量的角度 117 double get_angle( bool isVector = true) 118 { 119 angle = atan2( e.y - s.y, e.x - s.x); 120 if(!isVector && angle < 0) 121 angle += PI; 122 return angle; 123 } 124 125 //点在线段或直线上 1:点在直线上 2点在s,e所在矩形内 126 bool has(const Point &_Off, bool isSegment = false) const 127 { 128 bool ff = sgn( xmult( s, e, _Off), 0) == 0; 129 if( !isSegment) return ff; 130 return ff 131 && sgn(_Off.x - min(s.x, e.x), 0) >= 0 && sgn(_Off.x - max(s.x, e.x), 0) <= 0 132 && sgn(_Off.y - min(s.y, e.y), 0) >= 0 && sgn(_Off.y - max(s.y, e.y), 0) <= 0; 133 } 134 135 //点到直线/线段的距离 136 double dis(const Point &_Off, bool isSegment = false) 137 { 138 ///化为一般式 139 pton(); 140 //到直线垂足的距离 141 double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b); 142 //如果是线段判断垂足 143 if(isSegment) 144 { 145 double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b); 146 double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b); 147 double xb = max(s.x, e.x); 148 double yb = max(s.y, e.y); 149 double xs = s.x + e.x - xb; 150 double ys = s.y + e.y - yb; 151 if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps) 152 td = min( getdis(_Off,s), getdis(_Off,e)); 153 } 154 return fabs(td); 155 } 156 157 //关于直线对称的点 158 Point mirror(const Point &_Off) 159 { 160 ///注意先转为一般式 161 Point ret; 162 double d = a * a + b * b; 163 ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d; 164 ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d; 165 return ret; 166 } 167 //计算两点的中垂线 168 static Line ppline(const Point &_a,const Point &_b) 169 { 170 Line ret; 171 ret.s.x = (_a.x + _b.x) / 2; 172 ret.s.y = (_a.y + _b.y) / 2; 173 //一般式 174 ret.a = _b.x - _a.x; 175 ret.b = _b.y - _a.y; 176 ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x; 177 //两点式 178 if(fabs(ret.a) > eps) 179 { 180 ret.e.y = 0.0; 181 ret.e.x = - ret.c / ret.a; 182 if(ret.e == ret. s) 183 { 184 ret.e.y = 1e10; 185 ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a; 186 } 187 } 188 else 189 { 190 ret.e.x = 0.0; 191 ret.e.y = - ret.c / ret.b; 192 if(ret.e == ret. s) 193 { 194 ret.e.x = 1e10; 195 ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b; 196 } 197 } 198 return ret; 199 } 200 201 //------------直线和直线(向量)------------- 202 //向量向左边平移t的距离 203 Line& moveLine( double t) 204 { 205 Point of; 206 of = Point( -( e.y - s.y), e.x - s.x); 207 double dis = sqrt( of.x * of.x + of.y * of.y); 208 of.x= of.x * t / dis, of.y = of.y * t / dis; 209 s = s + of, e = e + of; 210 return *this; 211 } 212 //直线重合 213 static bool equal(const Line &_st,const Line &_se) 214 { 215 return _st.has( _se.e) && _se.has( _st.s); 216 } 217 //直线平行 218 static bool parallel(const Line &_st,const Line &_se) 219 { 220 return sgn( _st / _se, 0) == 0; 221 } 222 //两直线(线段)交点 223 //返回-1代表平行,0代表重合,1代表相交 224 static bool crossLPt(const Line &_st,const Line &_se, Point &ret) 225 { 226 if(parallel(_st,_se)) 227 { 228 if(Line::equal(_st,_se)) return 0; 229 return -1; 230 } 231 ret = _st.s; 232 double t = ( Line(_st.s,_se.s) / _se) / ( _st / _se); 233 ret.x += (_st.e.x - _st.s.x) * t; 234 ret.y += (_st.e.y - _st.s.y) * t; 235 return 1; 236 } 237 //------------线段和直线(向量)---------- 238 //直线和线段相交 239 //参数:直线[_st],线段[_se] 240 friend bool crossSL( Line &_st, Line &_se) 241 { 242 return sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.e), 0) >= 0; 243 } 244 245 //判断线段是否相交(注意添加eps) 246 static bool isCrossSS( const Line &_st, const Line &_se) 247 { 248 //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交 249 //2.跨立试验(等于0时端点重合) 250 return 251 max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) && 252 max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) && 253 max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) && 254 max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) && 255 sgn( xmult( _se.s, _st.s, _se.e) * xmult( _se.s, _se.e, _st.s), 0) >= 0 && 256 sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.s), 0) >= 0; 257 } 258 }; 259 260 //寻找凸包的graham 扫描法所需的排序函数 261 Point gsort; 262 bool gcmp( const Point &ta, const Point &tb)/// 选取与最后一条确定边夹角最小的点,即余弦值最大者 263 { 264 double tmp = xmult( gsort, ta, tb); 265 if( fabs( tmp) < eps) 266 return getdis( gsort, ta) < getdis( gsort, tb); 267 else if( tmp > 0) 268 return 1; 269 return 0; 270 } 271 272 class Polygon 273 { 274 public: 275 const static int maxpn = 5e4+7; 276 Point pt[maxpn];//点(顺时针或逆时针) 277 Line dq[maxpn]; //求半平面交打开注释 278 int n;//点的个数 279 280 281 //求多边形面积,多边形内点必须顺时针或逆时针 282 double area() 283 { 284 double ans = 0.0; 285 for(int i = 0; i < n; i ++) 286 { 287 int nt = (i + 1) % n; 288 ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y; 289 } 290 return fabs( ans / 2.0); 291 } 292 //求多边形重心,多边形内点必须顺时针或逆时针 293 Point gravity() 294 { 295 Point ans; 296 ans.x = ans.y = 0.0; 297 double area = 0.0; 298 for(int i = 0; i < n; i ++) 299 { 300 int nt = (i + 1) % n; 301 double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y; 302 area += tp; 303 ans.x += tp * (pt[i].x + pt[nt].x); 304 ans.y += tp * (pt[i].y + pt[nt].y); 305 } 306 ans.x /= 3 * area; 307 ans.y /= 3 * area; 308 return ans; 309 } 310 //判断点是否在任意多边形内[射线法],O(n) 311 bool ahas( Point &_Off) 312 { 313 int ret = 0; 314 double infv = 1e20;//坐标系最大范围 315 Line l = Line( _Off, Point( -infv ,_Off.y)); 316 for(int i = 0; i < n; i ++) 317 { 318 Line ln = Line( pt[i], pt[(i + 1) % n]); 319 if(fabs(ln.s.y - ln.e.y) > eps) 320 { 321 Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e; 322 if( ( fabs( tp.y - _Off.y) < eps && tp.x < _Off.x + eps) || Line::isCrossSS( ln, l)) 323 ret++; 324 } 325 else if( Line::isCrossSS( ln, l)) 326 ret++; 327 } 328 return ret&1; 329 } 330 331 //判断任意点是否在凸包内,O(logn) 332 bool bhas( Point & p) 333 { 334 if( n < 3) 335 return false; 336 if( xmult( pt[0], p, pt[1]) > eps) 337 return false; 338 if( xmult( pt[0], p, pt[n-1]) < -eps) 339 return false; 340 int l = 2,r = n-1; 341 int line = -1; 342 while( l <= r) 343 { 344 int mid = ( l + r) >> 1; 345 if( xmult( pt[0], p, pt[mid]) >= 0) 346 line = mid,r = mid - 1; 347 else l = mid + 1; 348 } 349 return xmult( pt[line-1], p, pt[line]) <= eps; 350 } 351 352 353 354 //凸多边形被直线分割 355 Polygon split( Line &_Off) 356 { 357 //注意确保多边形能被分割 358 Polygon ret; 359 Point spt[2]; 360 double tp = 0.0, np; 361 bool flag = true; 362 int i, pn = 0, spn = 0; 363 for(i = 0; i < n; i ++) 364 { 365 if(flag) 366 pt[pn ++] = pt[i]; 367 else 368 ret.pt[ret.n ++] = pt[i]; 369 np = xmult( _Off.s, _Off.e, pt[(i + 1) % n]); 370 if(tp * np < -eps) 371 { 372 flag = !flag; 373 Line::crossLPt( _Off, Line(pt[i], pt[(i + 1) % n]), spt[spn++]); 374 } 375 tp = (fabs(np) > eps)?np: tp; 376 } 377 ret.pt[ret.n ++] = spt[0]; 378 ret.pt[ret.n ++] = spt[1]; 379 n = pn; 380 return ret; 381 } 382 383 384 /** 卷包裹法求点集凸包,_p为输入点集,_n为点的数量 **/ 385 void ConvexClosure( Point _p[], int _n) 386 { 387 sort( _p, _p + _n); 388 n = 0; 389 for(int i = 0; i < _n; i++) 390 { 391 while( n > 1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0) 392 n--; 393 pt[n++] = _p[i]; 394 } 395 int _key = n; 396 for(int i = _n - 2; i >= 0; i--) 397 { 398 while( n > _key && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0) 399 n--; 400 pt[n++] = _p[i]; 401 } 402 if(n>1) n--;//除去重复的点,该点已是凸包凸包起点 403 } 404 /****** 寻找凸包的graham 扫描法********************/ 405 /****** _p为输入的点集,_n为点的数量****************/ 406 407 void graham( Point _p[], int _n) 408 { 409 int cur=0; 410 for(int i = 1; i < _n; i++) 411 if( sgn( _p[cur].y, _p[i].y) > 0 || ( sgn( _p[cur].y, _p[i].y) == 0 && sgn( _p[cur].x, _p[i].x) > 0) ) 412 cur = i; 413 swap( _p[cur], _p[0]); 414 n = 0, gsort = pt[n++] = _p[0]; 415 if( _n <= 1) return; 416 sort( _p + 1, _p+_n ,gcmp); 417 pt[n++] = _p[1]; 418 for(int i = 2; i < _n; i++) 419 { 420 while(n>1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)// 当凸包退化成直线时需特别注意n 421 n--; 422 pt[n++] = _p[i]; 423 } 424 } 425 //凸包旋转卡壳(注意点必须顺时针或逆时针排列) 426 //返回值凸包直径的平方(最远两点距离的平方) 427 pair<Point,Point> rotating_calipers() 428 { 429 int i = 1 % n; 430 double ret = 0.0; 431 pt[n] = pt[0]; 432 pair<Point,Point>ans=make_pair(pt[0],pt[0]); 433 for(int j = 0; j < n; j ++) 434 { 435 while( fabs( xmult( pt[i+1], pt[j], pt[j + 1])) > fabs( xmult( pt[i], pt[j], pt[j + 1])) + eps) 436 i = (i + 1) % n; 437 //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点 438 if(ret < getdis2(pt[i],pt[j])) ret = getdis2(pt[i],pt[j]), ans = make_pair(pt[i],pt[j]); 439 if(ret < getdis2(pt[i+1],pt[j+1])) ret = getdis(pt[i+1],pt[j+1]), ans = make_pair(pt[i+1],pt[j+1]); 440 } 441 return ans; 442 } 443 444 //凸包旋转卡壳(注意点必须逆时针排列) 445 //返回值两凸包的最短距离 446 double rotating_calipers( Polygon &_Off) 447 { 448 int i = 0; 449 double ret = 1e10;//inf 450 pt[n] = pt[0]; 451 _Off.pt[_Off.n] = _Off.pt[0]; 452 //注意凸包必须逆时针排列且pt[0]是左下角点的位置 453 while( _Off.pt[i + 1].y > _Off.pt[i].y) 454 i = (i + 1) % _Off.n; 455 for(int j = 0; j < n; j ++) 456 { 457 double tp; 458 //逆时针时为 >,顺时针则相反 459 while((tp = xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps) 460 i = (i + 1) % _Off.n; 461 //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段 462 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true)); 463 ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true)); 464 if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断 465 { 466 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true)); 467 ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true)); 468 } 469 } 470 return ret; 471 } 472 473 //-----------半平面交------------- 474 //复杂度:O(nlog2(n)) 475 //获取半平面交的多边形(多边形的核) 476 //参数:向量集合[l],向量数量[ln];(半平面方向在向量左边) 477 //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形(不存在区域或区域面积无穷大) 478 int judege( Line &_lx, Line &_ly, Line &_lz) 479 { 480 Point tmp; 481 Line::crossLPt(_lx,_ly,tmp); 482 return sgn(xmult(_lz.s,tmp,_lz.e),0); 483 } 484 int halfPanelCross(Line L[], int ln) 485 { 486 int i, tn, bot, top; 487 for(int i = 0; i < ln; i++) 488 L[i].get_angle(); 489 sort(L, L + ln); 490 //平面在向量左边的筛选 491 for(i = tn = 1; i < ln; i ++) 492 if(fabs(L[i].angle - L[i - 1].angle) > eps) 493 L[tn ++] = L[i]; 494 ln = tn, n = 0, bot = 0, top = 1; 495 dq[0] = L[0], dq[1] = L[1]; 496 for(i = 2; i < ln; i ++) 497 { 498 while(bot < top && judege(dq[top],dq[top-1],L[i]) > 0) 499 top --; 500 while(bot < top && judege(dq[bot],dq[bot+1],L[i]) > 0) 501 bot ++; 502 dq[++ top] = L[i]; 503 } 504 while(bot < top && judege(dq[top],dq[top-1],dq[bot]) > 0) 505 top --; 506 while(bot < top && judege(dq[bot],dq[bot+1],dq[top]) > 0) 507 bot ++; 508 //若半平面交退化为点或线 509 // if(top <= bot + 1) 510 // return 0; 511 dq[++top] = dq[bot]; 512 for(i = bot; i < top; i ++) 513 Line::crossLPt(dq[i],dq[i + 1],pt[n++]); 514 return n; 515 } 516 }; 517 518 Polygon pa,pb; 519 520 int main(void) 521 { 522 while(~scanf("%d%d",&pa.n,&pb.n)&&(pa.n||pb.n)) 523 { 524 for(int i=0;i<pa.n;i++) 525 scanf("%lf%lf",&pa.pt[i].x,&pa.pt[i].y); 526 for(int i=0;i<pb.n;i++) 527 scanf("%lf%lf",&pb.pt[i].x,&pb.pt[i].y); 528 printf("%.5f\n",pa.rotating_calipers(pb)); 529 } 530 return 0; 531 }
作者:weeping
出处:www.cnblogs.com/weeping/
本文版权归作者和博客园共有,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接,否则保留追究法律责任的权利。