poj3608 Bridge Across Islands

地址：http://poj.org/problem?id=3608

题目：

 

Bridge Across Islands
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.

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

POJ Founder Monthly Contest – 2008.06.29, Lei Tao

思路：

　　两凸包间的最近点对，套模板。

  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 }

 

posted @ 2017-10-12 13:44  weeping  阅读(369)  评论(0)  编辑  收藏  举报