地址：http://poj.org/problem?id=3525
题目：
Most Distant Point from the Sea
| Time Limit: 5000MS | Memory Limit: 65536K | |||
| Total Submissions: 5167 | Accepted: 2331 | Special Judge | ||
Description
The main land of Japan called Honshu is an island surrounded by the sea. In such an island, it is natural to ask a question: “Where is the most distant point from the sea?” The answer to this question for Honshu was found in 1996. The most distant point is located in former Usuda Town, Nagano Prefecture, whose distance from the sea is 114.86 km.
In this problem, you are asked to write a program which, given a map of an island, finds the most distant point from the sea in the island, and reports its distance from the sea. In order to simplify the problem, we only consider maps representable by convex polygons.
Input
The input consists of multiple datasets. Each dataset represents a map of an island, which is a convex polygon. The format of a dataset is as follows.
| n | ||
| x1 | y1 | |
| ⋮ | ||
| xn | yn |
Every input item in a dataset is a non-negative integer. Two input items in a line are separated by a space.
n in the first line is the number of vertices of the polygon, satisfying 3 ≤ n ≤ 100. Subsequent n lines are the x- and y-coordinates of the n vertices. Line segments (xi, yi)–(xi+1, yi+1) (1 ≤ i ≤ n − 1) and the line segment (xn, yn)–(x1, y1) form the border of the polygon in counterclockwise order. That is, these line segments see the inside of the polygon in the left of their directions. All coordinate values are between 0 and 10000, inclusive.
You can assume that the polygon is simple, that is, its border never crosses or touches itself. As stated above, the given polygon is always a convex one.
The last dataset is followed by a line containing a single zero.
Output
For each dataset in the input, one line containing the distance of the most distant point from the sea should be output. An output line should not contain extra characters such as spaces. The answer should not have an error greater than 0.00001 (10−5). You may output any number of digits after the decimal point, provided that the above accuracy condition is satisfied.
Sample Input
4
0 0
10000 0
10000 10000
0 10000
3
0 0
10000 0
7000 1000
6
0 40
100 20
250 40
250 70
100 90
0 70
3
0 0
10000 10000
5000 5001
0
Sample Output
5000.000000 494.233641 34.542948 0.353553
思路：二分距离d，然后将所有边内推距离d，之后用半平面交判断即可。
至于怎么内推，找出线段ab的方向向量（x,y）,那么向量c(-y,x)即是ab的垂线，再将c单位化后乘以d，线段(a+c)(b+c)即是内推后的线段位置。
1 #include <iostream>
2 #include <cstdio>
3 #include <cmath>
4 #include <algorithm>
5
6
7 using namespace std;
8 const double eps = 1e-10;
9 //点
10 class Point
11 {
12 public:
13 double x, y;
14
15 Point(){}
16 Point(double x, double y):x(x),y(y){}
17
18 bool operator < (const Point &_se) const
19 {
20 return x<_se.x || (x==_se.x && y<_se.y);
21 }
22 /*******判断ta与tb的大小关系*******/
23 static int sgn(double ta,double tb)
24 {
25 if(fabs(ta-tb)<eps)return 0;
26 if(ta<tb) return -1;
27 return 1;
28 }
29 static double xmult(const Point &po, const Point &ps, const Point &pe)
30 {
31 return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.y);
32 }
33 friend Point operator + (const Point &_st,const Point &_se)
34 {
35 return Point(_st.x + _se.x, _st.y + _se.y);
36 }
37 friend Point operator - (const Point &_st,const Point &_se)
38 {
39 return Point(_st.x - _se.x, _st.y - _se.y);
40 }
41 //点位置相同(double类型)
42 bool operator == (const Point &_off) const
43 {
44 return Point::sgn(x, _off.x) == 0 && Point::sgn(y, _off.y) == 0;
45 }
46 //点位置不同(double类型)
47 bool operator != (const Point &_Off) const
48 {
49 return ((*this) == _Off) == false;
50 }
51 //两点间距离的平方
52 static double dis2(const Point &_st,const Point &_se)
53 {
54 return (_st.x - _se.x) * (_st.x - _se.x) + (_st.y - _se.y) * (_st.y - _se.y);
55 }
56 //两点间距离
57 static double dis(const Point &_st, const Point &_se)
58 {
59 return sqrt((_st.x - _se.x) * (_st.x - _se.x) + (_st.y - _se.y) * (_st.y - _se.y));
60 }
61 };
62 //两点表示的向量
63 class Line
64 {
65 public:
66 Point s, e;//两点表示，起点[s]，终点[e]
67 double a, b, c;//一般式,ax+by+c=0
68
69 Line(){}
70 Line(const Point &s, const Point &e):s(s),e(e){}
71 Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){}
72
73 //向量与点的叉乘,参数：点[_Off]
74 //[点相对向量位置判断]
75 double operator /(const Point &_Off) const
76 {
77 return (_Off.y - s.y) * (e.x - s.x) - (_Off.x - s.x) * (e.y - s.y);
78 }
79 //向量与向量的叉乘,参数：向量[_Off]
80 friend double operator /(const Line &_st,const Line &_se)
81 {
82 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);
83 }
84 friend double operator *(const Line &_st,const Line &_se)
85 {
86 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);
87 }
88 //从两点表示转换为一般表示
89 //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2
90 bool pton()
91 {
92 a = e.y - s.y;
93 b = s.x - e.x;
94 c = e.x * s.y - e.y * s.x;
95 return true;
96 }
97
98 //-----------点和直线（向量）-----------
99 //点在向量左边（右边的小于号改成大于号即可,在对应直线上则加上=号）
100 //参数：点[_Off],向量[_Ori]
101 friend bool operator<(const Point &_Off, const Line &_Ori)
102 {
103 return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x)
104 < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x);
105 }
106
107 //点在直线上,参数：点[_Off]
108 bool lhas(const Point &_Off) const
109 {
110 return Point::sgn((*this) / _Off, 0) == 0;
111 }
112 //点在线段上,参数：点[_Off]
113 bool shas(const Point &_Off) const
114 {
115 return lhas(_Off)
116 && Point::sgn(_Off.x - min(s.x, e.x), 0) > 0 && Point::sgn(_Off.x - max(s.x, e.x), 0) < 0
117 && Point::sgn(_Off.y - min(s.y, e.y), 0) > 0 && Point::sgn(_Off.y - max(s.y, e.y), 0) < 0;
118 }
119
120 //点到直线/线段的距离
121 //参数： 点[_Off], 是否是线段[isSegment](默认为直线)
122 double dis(const Point &_Off, bool isSegment = false)
123 {
124 ///化为一般式
125 pton();
126
127 //到直线垂足的距离
128 double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b);
129
130 //如果是线段判断垂足
131 if(isSegment)
132 {
133 double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b);
134 double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b);
135 double xb = max(s.x, e.x);
136 double yb = max(s.y, e.y);
137 double xs = s.x + e.x - xb;
138 double ys = s.y + e.y - yb;
139 if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps)
140 td = min(Point::dis(_Off,s), Point::dis(_Off,e));
141 }
142
143 return fabs(td);
144 }
145
146 //关于直线对称的点
147 Point mirror(const Point &_Off) const
148 {
149 ///注意先转为一般式
150 Point ret;
151 double d = a * a + b * b;
152 ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d;
153 ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d;
154 return ret;
155 }
156 //计算两点的中垂线
157 static Line ppline(const Point &_a, const Point &_b)
158 {
159 Line ret;
160 ret.s.x = (_a.x + _b.x) / 2;
161 ret.s.y = (_a.y + _b.y) / 2;
162 //一般式
163 ret.a = _b.x - _a.x;
164 ret.b = _b.y - _a.y;
165 ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x;
166 //两点式
167 if(std::fabs(ret.a) > eps)
168 {
169 ret.e.y = 0.0;
170 ret.e.x = - ret.c / ret.a;
171 if(ret.e == ret. s)
172 {
173 ret.e.y = 1e10;
174 ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a;
175 }
176 }
177 else
178 {
179 ret.e.x = 0.0;
180 ret.e.y = - ret.c / ret.b;
181 if(ret.e == ret. s)
182 {
183 ret.e.x = 1e10;
184 ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b;
185 }
186 }
187 return ret;
188 }
189
190 //------------直线和直线（向量）-------------
191 //直线重合,参数：直线向量[_st],[_se]
192 static bool equal(const Line &_st, const Line &_se)
193 {
194 return _st.lhas(_se.e) && _se.lhas(_se.s);
195 }
196 //直线平行，参数：直线向量[_st],[_se]
197 static bool parallel(const Line &_st,const Line &_se)
198 {
199 return Point::sgn(_st / _se, 0) == 0;
200 }
201 //两直线（线段）交点，参数：直线向量[_st],[_se],交点
202 //返回-1代表平行，0代表重合，1代表相交
203 static bool crossLPt(const Line &_st,const Line &_se,Point &ret)
204 {
205 if(Line::parallel(_st,_se))
206 {
207 if(Line::equal(_st,_se)) return 0;
208 return -1;
209 }
210 ret = _st.s;
211 double t = (Line(_st.s,_se.s)/_se)/(_st/_se);
212 ret.x += (_st.e.x - _st.s.x) * t;
213 ret.y += (_st.e.y - _st.s.y) * t;
214 return 1;
215 }
216 //------------线段和直线（向量）----------
217 //线段和直线交
218 //参数：直线[_st],线段[_se]
219 friend bool crossSL(const Line &_st,const Line &_se)
220 {
221 return Point::sgn((_st / _se.s) * (_st / _se.e) ,0) <= 0;
222 }
223
224 //------------线段和线段（向量）----------
225 //判断线段是否相交(注意添加eps)，参数：线段[_st],线段[_se]
226 static bool isCrossSS(const Line &_st,const Line &_se)
227 {
228 //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
229 //2.跨立试验（等于0时端点重合）
230 return
231 max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) &&
232 max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) &&
233 max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) &&
234 max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) &&
235 Point::sgn((_st / Line(_st.s, _se.s)) * (_st / Line(_st.s, _se.e)), 0) <= 0 &&
236 Point::sgn((_se / Line(_se.s, _st.s)) * (_se / Line(_se.s, _st.e)), 0) <= 0;
237 }
238 };
239 class Polygon
240 {
241 public:
242 const static int maxpn = 100;
243 Point pt[maxpn];//点（顺时针或逆时针）
244 int n;//点的个数
245
246 Point& operator[](int _p)
247 {
248 return pt[_p];
249 }
250
251 //求多边形面积，多边形内点必须顺时针或逆时针
252 double area() const
253 {
254 double ans = 0.0;
255 for(int i = 0; i < n; i ++)
256 {
257 int nt = (i + 1) % n;
258 ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
259 }
260 return fabs(ans / 2.0);
261 }
262 //求多边形重心，多边形内点必须顺时针或逆时针
263 Point gravity() const
264 {
265 Point ans;
266 ans.x = ans.y = 0.0;
267 double area = 0.0;
268 for(int i = 0; i < n; i ++)
269 {
270 int nt = (i + 1) % n;
271 double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
272 area += tp;
273 ans.x += tp * (pt[i].x + pt[nt].x);
274 ans.y += tp * (pt[i].y + pt[nt].y);
275 }
276 ans.x /= 3 * area;
277 ans.y /= 3 * area;
278 return ans;
279 }
280 //判断点在凸多边形内，参数：点[_Off]
281 bool chas(const Point &_Off) const
282 {
283 double tp = 0, np;
284 for(int i = 0; i < n; i ++)
285 {
286 np = Line(pt[i], pt[(i + 1) % n]) / _Off;
287 if(tp * np < -eps)
288 return false;
289 tp = (fabs(np) > eps)?np: tp;
290 }
291 return true;
292 }
293 //判断点是否在任意多边形内[射线法]，O(n)
294 bool ahas(const Point &_Off) const
295 {
296 int ret = 0;
297 double infv = 1e-10;//坐标系最大范围
298 Line l = Line(_Off, Point( -infv ,_Off.y));
299 for(int i = 0; i < n; i ++)
300 {
301 Line ln = Line(pt[i], pt[(i + 1) % n]);
302 if(fabs(ln.s.y - ln.e.y) > eps)
303 {
304 Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e;
305 if(fabs(tp.y - _Off.y) < eps && tp.x < _Off.x + eps)
306 ret ++;
307 }
308 else if(Line::isCrossSS(ln,l))
309 ret ++;
310 }
311 return (ret % 2 == 1);
312 }
313 //凸多边形被直线分割,参数：直线[_Off]
314 Polygon split(Line _Off)
315 {
316 //注意确保多边形能被分割
317 Polygon ret;
318 Point spt[2];
319 double tp = 0.0, np;
320 bool flag = true;
321 int i, pn = 0, spn = 0;
322 for(i = 0; i < n; i ++)
323 {
324 if(flag)
325 pt[pn ++] = pt[i];
326 else
327 ret.pt[ret.n ++] = pt[i];
328 np = _Off / pt[(i + 1) % n];
329 if(tp * np < -eps)
330 {
331 flag = !flag;
332 Line::crossLPt(_Off,Line(pt[i], pt[(i + 1) % n]),spt[spn++]);
333 }
334 tp = (fabs(np) > eps)?np: tp;
335 }
336 ret.pt[ret.n ++] = spt[0];
337 ret.pt[ret.n ++] = spt[1];
338 n = pn;
339 return ret;
340 }
341
342
343 /** 卷包裹法求点集凸包，_p为输入点集，_n为点的数量 **/
344 void ConvexClosure(Point _p[],int _n)
345 {
346 sort(_p,_p+_n);
347 n=0;
348 for(int i=0;i<_n;i++)
349 {
350 while(n>1&&Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<=0)
351 n--;
352 pt[n++]=_p[i];
353 }
354 int _key=n;
355 for(int i=_n-2;i>=0;i--)
356 {
357 while(n>_key&&Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<=0)
358 n--;
359 pt[n++]=_p[i];
360 }
361 if(n>1) n--;//除去重复的点，该点已是凸包凸包起点
362 }
363 // /****** 寻找凸包的graham 扫描法********************/
364 // /****** _p为输入的点集,_n为点的数量****************/
365 // /**使用时需把gmp函数放在类外，并且看情况修改pt[0]**/
366 // bool gcmp(const Point &ta,const Point &tb)/// 选取与最后一条确定边夹角最小的点，即余弦值最大者
367 // {
368 // double tmp=Line(pt[0],ta)/Line(pt[0],tb);
369 // if(Point::sgn(tmp,0)==0)
370 // return Point::dis(pt[0],ta)<Point::dis(pt[0],tb);
371 // else if(tmp>0)
372 // return 1;
373 // return 0;
374 // }
375 // void graham(Point _p[],int _n)
376 // {
377 // int cur=0;
378 // for(int i=1;i<_n;i++)
379 // if(Point::sgn(_p[cur].y,_p[i].y)>0 || (Point::sgn(_p[cur].y,_p[i].y)==0 && Point::sgn(_p[cur].x,_p[i].x)>0))
380 // cur=i;
381 // swap(_p[cur],_p[0]);
382 // n=0,pt[n++]=_p[0];
383 // if(_n==1) return;
384 // sort(_p+1,_p+_n,Polygon::gcmp);
385 // pt[n++]=_p[1],pt[n++]=_p[2];
386 // for(int i=3;i<_n;i++)
387 // {
388 // while(Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<0)
389 // n--;
390 // pt[n++]=_p[i];
391 // }
392 // }
393 //凸包旋转卡壳(注意点必须顺时针或逆时针排列)
394 //返回值凸包直径的平方（最远两点距离的平方）
395 double rotating_calipers()
396 {
397 int i = 1;
398 double ret = 0.0;
399 pt[n] = pt[0];
400 for(int j = 0; j < n; j ++)
401 {
402 while(fabs(Point::xmult(pt[i+1],pt[j], pt[j + 1])) > fabs(Point::xmult(pt[i],pt[j], pt[j + 1])) + eps)
403 i = (i + 1) % n;
404 //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点
405 ret = (ret, max(Point::dis(pt[i],pt[j]), Point::dis(pt[i + 1],pt[j + 1])));
406 }
407 return ret;
408 }
409
410 //凸包旋转卡壳(注意点必须逆时针排列)
411 //返回值两凸包的最短距离
412 double rotating_calipers(Polygon &_Off)
413 {
414 int i = 0;
415 double ret = 1e10;//inf
416 pt[n] = pt[0];
417 _Off.pt[_Off.n] = _Off.pt[0];
418 //注意凸包必须逆时针排列且pt[0]是左下角点的位置
419 while(_Off.pt[i + 1].y > _Off.pt[i].y)
420 i = (i + 1) % _Off.n;
421 for(int j = 0; j < n; j ++)
422 {
423 double tp;
424 //逆时针时为 >,顺时针则相反
425 while((tp = Point::xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - Point::xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps)
426 i = (i + 1) % _Off.n;
427 //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段
428 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true));
429 ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true));
430 if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断
431 {
432 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true));
433 ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true));
434 }
435 }
436 return ret;
437 }
438
439 //-----------半平面交-------------
440 //复杂度:O(nlog2(n))
441 //#include <algorithm>
442 //半平面计算极角函数[如果考虑效率可以用成员变量记录]
443 static double hpc_pa(const Line &_Off)
444 {
445 return atan2(_Off.e.y - _Off.s.y, _Off.e.x - _Off.s.x);
446 }
447 //半平面交排序函数[优先顺序: 1.极角 2.前面的直线在后面的左边]
448 static bool hpc_cmp(const Line &l, const Line &r)
449 {
450 double lp = hpc_pa(l), rp = hpc_pa(r);
451 if(fabs(lp - rp) > eps)
452 return lp < rp;
453 return Point::xmult(r.s,l.s, r.e) < -eps;
454 }
455 static int judege(const Line &_lx,const Line &_ly,const Line &_lz)
456 {
457 Point tmp;
458 Line::crossLPt(_lx,_ly,tmp);
459 return Point::sgn(Point::xmult(_lz.s,tmp,_lz.e),0);
460 }
461 //获取半平面交的多边形（多边形的核）
462 //参数：向量集合[l]，向量数量[ln];(半平面方向在向量左边）
463 //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形（不存在区域或区域面积无穷大）
464 Polygon& halfPanelCross(Line _Off[], int ln)
465 {
466 Line dequeue[maxpn];//用于计算的双端队列
467 int i, tn, bot, top;
468 sort(_Off, _Off + ln, hpc_cmp);
469 //平面在向量左边的筛选
470 for(i = tn = 1; i < ln; i ++)
471 if(fabs(hpc_pa(_Off[i]) - hpc_pa(_Off[i - 1])) > eps)
472 _Off[tn ++] = _Off[i];
473 ln = tn, n = 0, bot = 0, top = 1;
474 dequeue[0] = _Off[0];
475 dequeue[1] = _Off[1];
476 for(i = 2; i < ln; i ++)
477 {
478 while(bot < top && Polygon::judege(dequeue[top],dequeue[top-1],_Off[i]) > 0)
479 top --;
480 while(bot < top && Polygon::judege(dequeue[bot],dequeue[bot+1],_Off[i]) > 0)
481 bot ++;
482 dequeue[++ top] = _Off[i];
483 }
484 while(bot < top && Polygon::judege(dequeue[top],dequeue[top-1],dequeue[bot]) > 0)
485 top --;
486 while(bot < top && Polygon::judege(dequeue[bot],dequeue[bot+1],dequeue[top]) > 0)
487 bot ++;
488 //计算交点(注意不同直线形成的交点可能重合)
489 if(top <= bot + 1)
490 return (*this);
491 for(i = bot; i < top; i ++)
492 Line::crossLPt(dequeue[i],dequeue[i + 1],pt[n++]);
493 if(bot < top + 1)
494 Line::crossLPt(dequeue[bot],dequeue[top],pt[n++]);
495 return (*this);
496 }
497 };
498
499 Point pt[100];
500 Line ln[100],ls[100];
501 Polygon py;
502 int n;
503
504 void ml(Line &lx,Line &ly,double t)
505 {
506 Point tmp=lx.e-lx.s,of;
507 of=Point(-tmp.y,tmp.x);
508 double dis=sqrt(of.x*of.x+of.y*of.y);
509 of.x=of.x*t/dis,of.y=of.y*t/dis;
510 ly.s=lx.s+of,ly.e=lx.e+of;
511 }
512 int sc(double x)
513 {
514 for(int i=0;i<n;i++)
515 ml(ln[i],ls[i],x);
516 //for(int i=0;i<n;i++)
517 // printf("%.2f %.2f %.2f %.2f
",ls[i].s.x,ls[i].s.y,ls[i].e.x,ls[i].e.y);
518 py.halfPanelCross(ls,n);
519 //printf("==%d
",py.n);
520 //for(int i=0;i<py.n;i++)
521 // printf("%.2f %.2f
",py.pt[i].x,py.pt[i].y);
522 return py.n;
523 }
524 int main(void)
525 {
526 while(scanf("%d",&n)&&n)
527 {
528 for(int i=0;i<n;i++)
529 scanf("%lf%lf",&pt[i].x,&pt[i].y);
530 pt[n]=pt[0];
531 for(int i=0;i<n;i++)
532 ln[i]=Line(pt[i],pt[i+1]);
533 //for(int i=0;i<n;i++)
534 // printf("%.2f %.2f %.2f %.2f
",ln[i].s.x,ln[i].s.y,ln[i].e.x,ln[i].e.y);
535 double l=0,r=1e8;
536 while(fabs(r-l)>=1e-5)
537 {
538 double mid=(r+l)/2.0;
539 if(sc(mid)) l=mid;
540 else r=mid;
541 }
542 printf("%.6f
",r);
543 }
544 return 0;
545 }