POJ3525 半平面交

题意：求某凸多边形内部离边界最远的点到边界的距离

首先介绍半平面、半平面交的概念：

半平面：对于一条有向直线，它的方向的左手侧就是它所划定的半平面范围。如图所示：

半平面交：多个半平面的交集。有点类似二元函数的线性规划。如图

求半平面交：用的kuangbin模板= =

sol：二分答案

　　二分距离值，按这个值把边界向内缩，再求新的半平面交。如图：

绿色的是原图形，蓝色是按距离值向里面缩进去之后得到的新图形。对这个新图做半平面交即可。

若半平面交存在，说明与边界的距离是该值的点存在（半平面交里面的点都是）

那么如何把图形向里面缩呢？其实在纸上画画就容易多了：

因为半平面是一条有向直线，所以只要有直线上的任意一点和方向就可以了，这个点在不在对应的那条边上都无所谓。

如图，棕色是一开始的边，黑色是它对应的半平面。

设当前二分答案值是mid，对这条边做一个长度为mid的法向量（图中绿色）

棕色边起点+法向量即得到新半平面上的一个点（图中蓝色）

新半平面的方向和原来一样，直接平移过来就行了。

在做法向量的时候采用向量的概念能简化问题。这里在kuangbin模板的基础上加上了向量类Vector：

 1 struct Vector:public point
 2 {
 3     Vector(){}
 4     Vector(double a,double b)
 5     {
 6         x=a;    y=b;
 7     }
 8     Vector(point _a,point _b)   //a->b
 9     {
10         double dx=_b.x-_a.x;
11         double dy=_b.y-_a.y;
12         x=dx;   y=dy;
13     }
14     Vector(line v)
15     {
16         double dx=v.b.x-v.a.x;
17         double dy=v.b.y-v.a.y;
18         x=dx;   y=dy;
19     }
20     double length()
21     {
22         return (sqrt(x*x+y*y));
23     }
24     Vector Normal()      //返回this的单位长度的法向量
25     {
26         double L=sqrt(x*x+y*y);
27         Vector Vans=Vector(-y/L,x/L);
28         return Vans;
29     }
30 };
31

重载了三种构造函数：给出向量的坐标值(x,y)、给出向量首尾两点、用一条线段作为向量

----------------------------------------------

AC Code：

  1 #include<vector>
  2 #include<list>
  3 #include<map>
  4 #include<set>
  5 #include<deque>
  6 #include<queue>
  7 #include<stack>
  8 #include<bitset>
  9 #include<algorithm>
 10 #include<functional>
 11 #include<numeric>
 12 #include<utility>
 13 #include<iostream>
 14 #include<sstream>
 15 #include<iomanip>
 16 #include<cstdio>
 17 #include<cmath>
 18 #include<cstdlib>
 19 #include<cctype>
 20 #include<string>
 21 #include<cstring>
 22 #include<cstdio>
 23 #include<cmath>
 24 #include<cstdlib>
 25 #include<ctime>
 26 #include<climits>
 27 #include<complex>
 28 #define mp make_pair
 29 #define pb push_back
 30 using namespace std;
 31 const double eps=1e-6;
 32 const double pi=acos(-1.0);
 33 const double inf=1e20;
 34 const int maxp=1111;
 35 int dblcmp(double d)
 36 {
 37     if (fabs(d)<eps)return 0;
 38     return d>eps?1:-1;
 39 }
 40 inline double sqr(double x){return x*x;}
 41 struct point
 42 {
 43     double x,y;
 44     point(){}
 45     point(double _x,double _y):
 46     x(_x),y(_y){};
 47     void input()
 48     {
 49         scanf("%lf%lf",&x,&y);
 50     }
 51     void output()
 52     {
 53         printf("%.2f %.2f\n",x,y);
 54     }
 55     bool operator==(point a)const
 56     {
 57         return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0;
 58     }
 59     bool operator<(point a)const
 60     {
 61         return dblcmp(a.x-x)==0?dblcmp(y-a.y)<0:x<a.x;
 62     }
 63     double len()
 64     {
 65         return hypot(x,y);
 66     }
 67     double len2()
 68     {
 69         return x*x+y*y;
 70     }
 71     double distance(point p)
 72     {
 73         return hypot(x-p.x,y-p.y);
 74     }
 75     point add(point p)
 76     {
 77         return point(x+p.x,y+p.y);
 78     }
 79     point sub(point p)
 80     {
 81         return point(x-p.x,y-p.y);
 82     }
 83     point mul(double b)
 84     {
 85         return point(x*b,y*b);
 86     }
 87     point div(double b)
 88     {
 89         return point(x/b,y/b);
 90     }
 91     double dot(point p)
 92     {
 93         return x*p.x+y*p.y;
 94     }
 95     double det(point p)
 96     {
 97         return x*p.y-y*p.x;
 98     }
 99     double rad(point a,point b)
100     {
101         point p=*this;
102         return fabs(atan2(fabs(a.sub(p).det(b.sub(p))),a.sub(p).dot(b.sub(p))));
103     }
104     point trunc(double r)
105     {
106         double l=len();
107         if (!dblcmp(l))return *this;
108         r/=l;
109         return point(x*r,y*r);
110     }
111     point rotleft()
112     {
113         return point(-y,x);
114     }
115     point rotright()
116     {
117         return point(y,-x);
118     }
119     point rotate(point p,double angle)//绕点p逆时针旋转angle角度
120     {
121         point v=this->sub(p);
122         double c=cos(angle),s=sin(angle);
123         return point(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);
124     }
125 };
126 struct line
127 {
128     point a,b;
129     line(){}
130     line(point _a,point _b)
131     {
132         a=_a;
133         b=_b;
134     }
135     bool operator==(line v)
136     {
137         return (a==v.a)&&(b==v.b);
138     }
139     //倾斜角angle
140     line(point p,double angle)
141     {
142         a=p;
143         if (dblcmp(angle-pi/2)==0)
144         {
145             b=a.add(point(0,1));
146         }
147         else
148         {
149             b=a.add(point(1,tan(angle)));
150         }
151     }
152     //ax+by+c=0
153     line(double _a,double _b,double _c)
154     {
155         if (dblcmp(_a)==0)
156         {
157             a=point(0,-_c/_b);
158             b=point(1,-_c/_b);
159         }
160         else if (dblcmp(_b)==0)
161         {
162             a=point(-_c/_a,0);
163             b=point(-_c/_a,1);
164         }
165         else
166         {
167             a=point(0,-_c/_b);
168             b=point(1,(-_c-_a)/_b);
169         }
170     }
171     void input()
172     {
173         a.input();
174         b.input();
175     }
176     void adjust()
177     {
178         if (b<a)swap(a,b);
179     }
180     double length()
181     {
182         return a.distance(b);
183     }
184     double angle()//直线倾斜角 0<=angle<180
185     {
186         double k=atan2(b.y-a.y,b.x-a.x);
187         if (dblcmp(k)<0)k+=pi;
188         if (dblcmp(k-pi)==0)k-=pi;
189         return k;
190     }
191     //点和线段关系
192     //1 在逆时针
193     //2 在顺时针
194     //3 平行
195     int relation(point p)
196     {
197         int c=dblcmp(p.sub(a).det(b.sub(a)));
198         if (c<0)return 1;
199         if (c>0)return 2;
200         return 3;
201     }
202     bool pointonseg(point p)
203     {
204         return dblcmp(p.sub(a).det(b.sub(a)))==0&&dblcmp(p.sub(a).dot(p.sub(b)))<=0;
205     }
206     bool parallel(line v)
207     {
208         return dblcmp(b.sub(a).det(v.b.sub(v.a)))==0;
209     }
210     //2 规范相交
211     //1 非规范相交
212     //0 不相交
213     int segcrossseg(line v)
214     {
215         int d1=dblcmp(b.sub(a).det(v.a.sub(a)));
216         int d2=dblcmp(b.sub(a).det(v.b.sub(a)));
217         int d3=dblcmp(v.b.sub(v.a).det(a.sub(v.a)));
218         int d4=dblcmp(v.b.sub(v.a).det(b.sub(v.a)));
219         if ((d1^d2)==-2&&(d3^d4)==-2)return 2;
220         return (d1==0&&dblcmp(v.a.sub(a).dot(v.a.sub(b)))<=0||
221                 d2==0&&dblcmp(v.b.sub(a).dot(v.b.sub(b)))<=0||
222                 d3==0&&dblcmp(a.sub(v.a).dot(a.sub(v.b)))<=0||
223                 d4==0&&dblcmp(b.sub(v.a).dot(b.sub(v.b)))<=0);
224     }
225     int linecrossseg(line v)//*this seg v line
226     {
227         int d1=dblcmp(b.sub(a).det(v.a.sub(a)));
228         int d2=dblcmp(b.sub(a).det(v.b.sub(a)));
229         if ((d1^d2)==-2)return 2;
230         return (d1==0||d2==0);
231     }
232     //0 平行
233     //1 重合
234     //2 相交
235     int linecrossline(line v)
236     {
237         if ((*this).parallel(v))
238         {
239             return v.relation(a)==3;
240         }
241         return 2;
242     }
243     point crosspoint(line v)
244     {
245         double a1=v.b.sub(v.a).det(a.sub(v.a));
246         double a2=v.b.sub(v.a).det(b.sub(v.a));
247         return point((a.x*a2-b.x*a1)/(a2-a1),(a.y*a2-b.y*a1)/(a2-a1));
248     }
249     double dispointtoline(point p)
250     {
251         return fabs(p.sub(a).det(b.sub(a)))/length();
252     }
253     double dispointtoseg(point p)
254     {
255         if (dblcmp(p.sub(b).dot(a.sub(b)))<0||dblcmp(p.sub(a).dot(b.sub(a)))<0)
256         {
257             return min(p.distance(a),p.distance(b));
258         }
259         return dispointtoline(p);
260     }
261     point lineprog(point p)
262     {
263         return a.add(b.sub(a).mul(b.sub(a).dot(p.sub(a))/b.sub(a).len2()));
264     }
265     point symmetrypoint(point p)
266     {
267         point q=lineprog(p);
268         return point(2*q.x-p.x,2*q.y-p.y);
269     }
270 };
271 
272 struct Vector:public point
273 {
274     Vector(){}
275     Vector(double a,double b)
276     {
277         x=a;    y=b;
278     }
279     Vector(point _a,point _b)   //a->b
280     {
281         double dx=_b.x-_a.x;
282         double dy=_b.y-_a.y;
283         x=dx;   y=dy;
284     }
285     Vector(line v)
286     {
287         double dx=v.b.x-v.a.x;
288         double dy=v.b.y-v.a.y;
289         x=dx;   y=dy;
290     }
291     double length()
292     {
293         return (sqrt(x*x+y*y));
294     }
295     Vector Normal()
296     {
297         double L=sqrt(x*x+y*y);
298         Vector Vans=Vector(-y/L,x/L);
299         return Vans;
300     }
301 };
302 
303 struct halfplane:public line    //半平面
304 {
305     double angle;
306     halfplane(){}
307     //表示向量 a->b逆时针(左侧)的半平面
308     halfplane(point _a,point _b)
309     {
310         a=_a;
311         b=_b;
312     }
313     halfplane(line v)
314     {
315         a=v.a;
316         b=v.b;
317     }
318     void calcangle()
319     {
320         angle=atan2(b.y-a.y,b.x-a.x);
321     }
322     bool operator<(const halfplane &b)const
323     {
324         return angle<b.angle;
325     }
326 };
327 struct halfplanes   //半平面交
328 {
329     int n;
330     halfplane hp[maxp];
331     point p[maxp];
332     int que[maxp];
333     int st,ed;
334     void push(halfplane tmp)
335     {
336         hp[n++]=tmp;
337     }
338     void unique()
339     {
340         int m=1,i;
341         for (i=1;i<n;i++)
342         {
343             if (dblcmp(hp[i].angle-hp[i-1].angle))hp[m++]=hp[i];
344             else if (dblcmp(hp[m-1].b.sub(hp[m-1].a).det(hp[i].a.sub(hp[m-1].a))>0))hp[m-1]=hp[i];
345         }
346         n=m;
347     }
348     bool halfplaneinsert()
349     {
350         int i;
351         for (i=0;i<n;i++)hp[i].calcangle();
352         sort(hp,hp+n);
353         unique();
354         que[st=0]=0;
355         que[ed=1]=1;
356         p[1]=hp[0].crosspoint(hp[1]);
357         for (i=2;i<n;i++)
358         {
359             while (st<ed&&dblcmp((hp[i].b.sub(hp[i].a).det(p[ed].sub(hp[i].a))))<0)ed--;
360             while (st<ed&&dblcmp((hp[i].b.sub(hp[i].a).det(p[st+1].sub(hp[i].a))))<0)st++;
361             que[++ed]=i;
362             if (hp[i].parallel(hp[que[ed-1]]))return false;
363             p[ed]=hp[i].crosspoint(hp[que[ed-1]]);
364         }
365         while (st<ed&&dblcmp(hp[que[st]].b.sub(hp[que[st]].a).det(p[ed].sub(hp[que[st]].a)))<0)ed--;
366         while (st<ed&&dblcmp(hp[que[ed]].b.sub(hp[que[ed]].a).det(p[st+1].sub(hp[que[ed]].a)))<0)st++;
367         if (st+1>=ed)return false;
368         return true;
369     }
370     /*
371     void getconvex(polygon &con)
372     {
373         p[st]=hp[que[st]].crosspoint(hp[que[ed]]);
374         con.n=ed-st+1;
375         int j=st,i=0;
376         for (;j<=ed;i++,j++)
377         {
378             con.p[i]=p[j];
379         }
380     }*/
381 };
382 
383 point p[1000];
384 Vector V[1000],Vt[1000];
385 halfplanes TH;
386 int n;
387 
388 int main()
389 {
390     //freopen("in.txt","r",stdin);
391 
392     while (cin>>n)
393     {
394         if (n==0)   break;
395         for (int i=0;i<n;i++)   //n points:[0..n-1]
396             p[i].input();
397 
398         for (int i=0;i<n;i++)   //v[i]:p[i]->p[i+1]
399         {
400             V[i]=Vector(p[i],p[(i+1)%n]);
401             //printf("vector:   %.6lf  %.6lf\n",V[i].x,V[i].y);
402             Vt[i]=V[i].Normal();
403         }
404 
405         double l=0,r=5001;
406         while (r-l>eps)
407         {
408             double mid=(l+r)/2;
409             //double mid=l+(r-l)/2;
410             TH.n=0;
411             //halfplanes TH;
412 
413             for (int i=0;i<n;i++)
414             {
415                 point t1=p[i].add(Vt[i].mul(mid));
416                 point t2=t1.add(V[i]);
417 
418                 line tmp=line(t1,t2);
419                 TH.push(halfplane(tmp));
420             }
421 
422             //printf("%.6lf  %d",mid,TH.n);
423 
424             if (TH.halfplaneinsert())
425             {
426                 //cout<<"OK"<<endl;
427                 l=mid;      //l=mid+0.00001;
428             }
429             else
430             {
431                 //cout<<"NO"<<endl;
432                 r=mid;      //r=mid-0.00001;
433             }
434         }
435         printf("%.6lf\n",l);
436     }
437     return 0;
438 }

View Code

推荐iPad上一个免费又好用的绘图板：Sketches

posted on 2015-02-05 20:47 Pentium.Labs 阅读(643) 评论(0) 编辑收藏举报

刷新页面返回顶部

Pentium.Labs

导航

公告

POJ3525 半平面交