Art Gallery

Art Gallery

http://poj.org/problem?id=1279

Time Limit: 1000MS   Memory Limit: 10000K
Total Submissions: 9684   Accepted: 3747

Description

The art galleries of the new and very futuristic building of the Center for Balkan Cooperation have the form of polygons (not necessarily convex). When a big exhibition is organized, watching over all of the pictures is a big security concern. Your task is that for a given gallery to write a program which finds the surface of the area of the floor, from which each point on the walls of the gallery is visible. On the figure 1. a map of a gallery is given in some co-ordinate system. The area wanted is shaded on the figure 2. 

Input

The number of tasks T that your program have to solve will be on the first row of the input file. Input data for each task start with an integer N, 5 <= N <= 1500. Each of the next N rows of the input will contain the co-ordinates of a vertex of the polygon ? two integers that fit in 16-bit integer type, separated by a single space. Following the row with the co-ordinates of the last vertex for the task comes the line with the number of vertices for the next test and so on.

Output

For each test you must write on one line the required surface - a number with exactly two digits after the decimal point (the number should be rounded to the second digit after the decimal point).

Sample Input

1
7
0 0
4 4
4 7
9 7
13 -1
8 -6
4 -4

Sample Output

80.00

 

 

要特判面积不存在的情况

  1 #include<iostream>
  2 #include<cstring>
  3 #include<cstdio>
  4 #include<vector>
  5 #include<cmath>
  6 #include<algorithm>
  7 using namespace std;
  8 const double eps=1e-8;
  9 const double INF=1e20;
 10 const double PI=acos(-1.0);
 11 
 12 
 13 int sgn(double x){
 14     if(fabs(x)<eps) return 0;
 15     if(x<0) return -1;
 16     else return 1;
 17 }
 18 inline double sqr(double x){return x*x;}
 19 struct Point{
 20     double x,y;
 21     Point(){}
 22     Point(double _x,double _y){
 23         x=_x;
 24         y=_y;
 25     }
 26     void input(){
 27         scanf("%lf %lf",&x,&y);
 28     }
 29     void output(){
 30         printf("%.2f %.2f\n",x,y);
 31     }
 32     bool operator == (const Point &b)const{
 33         return sgn(x-b.x) == 0 && sgn(y-b.y)== 0;
 34     }
 35     bool operator < (const Point &b)const{
 36         return sgn(x-b.x)==0?sgn(y-b.y)<0:x<b.x;
 37     }
 38     Point operator - (const Point &b)const{
 39         return Point(x-b.x,y-b.y);
 40     }
 41     //叉积
 42     double operator ^ (const Point &b)const{
 43         return x*b.y-y*b.x;
 44     }
 45     //点积
 46     double operator * (const Point &b)const{
 47         return x*b.x+y*b.y;
 48     }
 49     //返回长度
 50     double len(){
 51         return hypot(x,y);
 52     }
 53     //返回长度的平方
 54     double len2(){
 55         return x*x+y*y;
 56     }
 57     //返回两点的距离
 58     double distance(Point p){
 59         return hypot(x-p.x,y-p.y);
 60     }
 61     Point operator + (const Point &b)const{
 62         return Point(x+b.x,y+b.y);
 63     }
 64     Point operator * (const double &k)const{
 65         return Point(x*k,y*k);
 66     }
 67     Point operator / (const double &k)const{
 68         return Point(x/k,y/k);
 69     }
 70 
 71     //计算pa和pb的夹角
 72     //就是求这个点看a,b所成的夹角
 73     ///LightOJ1202
 74     double rad(Point a,Point b){
 75         Point p=*this;
 76         return fabs(atan2(fabs((a-p)^(b-p)),(a-p)*(b-p)));
 77     }
 78     //化为长度为r的向量
 79     Point trunc(double r){
 80         double l=len();
 81         if(!sgn(l)) return *this;
 82         r/=l;
 83         return Point(x*r,y*r);
 84     }
 85     //逆时针转90度
 86     Point rotleft(){
 87         return Point(-y,x);
 88     }
 89     //顺时针转90度
 90     Point rotright(){
 91         return Point(y,-x);
 92     }
 93     //绕着p点逆时针旋转angle
 94     Point rotate(Point p,double angle){
 95         Point v=(*this) -p;
 96         double c=cos(angle),s=sin(angle);
 97         return Point(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);
 98     }
 99 };
100 
101 struct Line{
102     Point s,e;
103     Line(){}
104     Line(Point _s,Point _e){
105         s=_s;
106         e=_e;
107     }
108     bool operator==(Line v){
109         return (s==v.s)&&(e==v.e);
110     }
111     //根据一个点和倾斜角angle确定直线，0<=angle<pi
112     Line(Point p,double angle){
113         s=p;
114         if(sgn(angle-PI/2)==0){
115             e=(s+Point(0,1));
116         }
117         else{
118             e=(s+Point(1,tan(angle)));
119         }
120     }
121     //ax+by+c=0;
122     Line(double a,double b,double c){
123         if(sgn(a)==0){
124             s=Point(0,-c/b);
125             e=Point(1,-c/b);
126         }
127         else if(sgn(b)==0){
128             s=Point(-c/a,0);
129             e=Point(-c/a,1);
130         }
131         else{
132             s=Point(0,-c/b);
133             e=Point(1,(-c-a)/b);
134         }
135     }
136     void input(){
137         s.input();
138         e.input();
139     }
140     void adjust(){
141         if(e<s) swap(s,e);
142     }
143     //求线段长度
144     double length(){
145         return s.distance(e);
146     }
147     //返回直线倾斜角 0<=angle<pi
148     double angle(){
149         double k=atan2(e.y-s.y,e.x-s.x);
150         if(sgn(k)<0) k+=PI;
151         if(sgn(k-PI)==0) k-=PI;
152         return k;
153     }
154     //点和直线的关系
155     //1 在左侧
156     //2 在右侧
157     //3 在直线上
158     int relation(Point p){
159         int c=sgn((p-s)^(e-s));
160         if(c<0) return 1;
161         else if(c>0) return 2;
162         else return 3;
163     }
164     //点在线段上的判断
165     bool pointonseg(Point p){
166         return sgn((p-s)^(e-s))==0&&sgn((p-s)*(p-e))<=0;
167     }
168     //两向量平行(对应直线平行或重合)
169     bool parallel(Line v){
170         return sgn((e-s)^(v.e-v.s))==0;
171     }
172     //两线段相交判断
173     //2 规范相交
174     //1 非规范相交
175     //0 不相交
176     int segcrossseg(Line v){
177         int d1=sgn((e-s)^(v.s-s));
178         int d2=sgn((e-s)^(v.e-s));
179         int d3=sgn((v.e-v.s)^(s-v.s));
180         int d4=sgn((v.e-v.s)^(e-v.s));
181         if((d1^d2)==-2&&(d3^d4)==-2) return 2;
182         return (d1==0&&sgn((v.s-s)*(v.s-e))<=0||
183                 d2==0&&sgn((v.e-s)*(v.e-e))<=0||
184                 d3==0&&sgn((s-v.s)*(s-v.e))<=0||
185                 d4==0&&sgn((e-v.s)*(e-v.e))<=0);
186     }
187     //直线和线段相交判断
188     //-*this line -v seg
189     //2 规范相交
190     //1 非规范相交
191     //0 不相交
192     int linecrossseg(Line v){
193         int d1=sgn((e-s)^(v.s-s));
194         int d2=sgn((e-s)^(v.e-s));
195         if((d1^d2)==-2) return 2;
196         return (d1==0||d2==0);
197     }
198     //两直线关系
199     //0 平行
200     //1 重合
201     //2 相交
202     int linecrossline(Line v){
203         if((*this).parallel(v))
204             return v.relation(s)==3;
205         return 2;
206     }
207     //求两直线的交点
208     //要保证两直线不平行或重合
209     Point crosspoint(Line v){
210         double a1=(v.e-v.s)^(s-v.s);
211         double a2=(v.e-v.s)^(e-v.s);
212         return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
213     }
214     //点到直线的距离
215     double dispointtoline(Point p){
216         return fabs((p-s)^(e-s))/length();
217     }
218     //点到线段的距离
219     double dispointtoseg(Point p){
220         if(sgn((p-s)*(e-s))<0||sgn((p-e)*(s-e))<0)
221             return min(p.distance(s),p.distance(e));
222         return dispointtoline(p);
223     }
224     //返回线段到线段的距离
225     //前提是两线段不相交，相交距离就是0了
226     double dissegtoseg(Line v){
227         return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));
228     }
229     //返回点P在直线上的投影
230     Point lineprog(Point p){
231         return s+(((e-s)*((e-s)*(p-s)))/((e-s).len2()));
232     }
233     //返回点P关于直线的对称点
234     Point symmetrypoint(Point p){
235         Point q=lineprog(p);
236         return Point(2*q.x-p.x,2*q.y-p.y);
237     }
238 };
239 
240 struct circle{
241     Point p;
242     double r;
243     circle(){}
244     circle(Point _p,double _r){
245         p=_p;
246         r=_r;
247     }
248 
249     circle(double x,double y,double _r){
250         p=Point(x,y);
251         r=_r;
252     }
253 
254     circle(Point a,Point b,Point c){///三角形的外接圆
255         Line u=Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
256         Line v=Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
257         p=u.crosspoint(v);
258         r=p.distance(a);
259     }
260 
261     circle(Point a,Point b,Point c,bool t){///三角形的内切圆
262         Line u,v;
263         double m=atan2(b.y-a.y,b.x-a.x),n=atan2(c.y-a.y,c.x-a.x);
264         u.s=a;
265         u.e=u.s+Point(cos((n+m)/2),sin((n+m)/2));
266         v.s=b;
267         m=atan2(a.y-b.y,a.x-b.x),n=atan2(c.y-b.y,c.x-b.x);
268         v.e=v.s+Point(cos((n+m)/2),sin((n+m)/2));
269         p=u.crosspoint(v);
270         r=Line(a,b).dispointtoseg(p);
271     }
272 
273     void input(){
274         p.input();
275         scanf("%lf",&r);
276     }
277 
278     void output(){
279         printf("%.2f %.2f %.2f\n",p.x,p.y,r);
280     }
281 
282     bool operator==(circle v){
283         return (p==v.p)&&sgn(r-v.r)==0;
284     }
285 
286     bool operator<(circle v)const{
287         return ((p<v.p)||((p==v.p)&&sgn(r-v.r)<0));
288     }
289 
290     double area(){
291         return PI*r*r;
292     }
293 
294     double circumference(){ ///周长
295         return 2*PI*r;
296     }
297 
298     int relation(Point b){///点和圆的关系  0圆外  1圆上  2圆内
299         double dst=b.distance(p);
300         if(sgn(dst-r)<0) return 2;
301         else if(sgn(dst-r)==0) return 1;
302         return 0;
303     }
304 
305     int relationseg(Line v){///线段和圆的关系，比较的是圆心到线段的距离和半径的关系
306         double dst=v.dispointtoseg(p);
307         if(sgn(dst-r)<0) return 2;
308         else if(sgn(dst-r)==0) return 1;
309         return 0;
310     }
311 
312     int relationline(Line v){///直线和圆的关系，比较的是圆心到直线的距离和半径的关系
313         double dst=v.dispointtoline(p);
314         if(sgn(dst-r)<0) return 2;
315         else if(sgn(dst-r)==0) return 1;
316         return 0;
317     }
318 
319     int relationcircle(circle v){///两圆的关系  5相离 4外切 3相交 2内切 1内含
320         double d=p.distance(v.p);
321         if(sgn(d-r-v.r)>0) return 5;
322         if(sgn(d-r-v.r)==0) return 4;
323         double l=fabs(r-v.r);
324         if(sgn(d-r-v.r)<0&&sgn(d-l)>0) return 3;
325         if(sgn(d-l)==0) return 2;
326         if(sgn(d-l)<0)  return 1;
327     }
328 
329     int pointcrosscircle(circle v,Point &p1,Point &p2){///求两个圆的交点，0没有交点 1一个交点 2两个交点
330         int rel=relationcircle(v);
331         if(rel == 1 || rel == 5) return 0;
332         double d=p.distance(v.p);
333         double l=(d*d+r*r-v.r*v.r)/2*d;
334         double h=sqrt(r*r-l*l);
335         Point tmp=p+(v.p-p).trunc(l);
336         p1=tmp+((v.p-p).rotleft().trunc(h));
337         p2=tmp+((v.p-p).rotright().trunc(h));
338         if(rel == 2 || rel == 4) return 1;
339         return 2;
340     }
341 
342     int pointcrossline(Line v,Point &p1,Point &p2){///求直线和圆的交点，返回交点的个数
343         if(!(*this).relationline(v)) return 0;
344         Point a=v.lineprog(p);
345         double d=v.dispointtoline(p);
346         d=sqrt(r*r-d*d);
347         if(sgn(d)==0) {
348             p1=a;
349             p2=a;
350             return 1;
351         }
352         p1=a+(v.e-v.s).trunc(d);
353         p2=a-(v.e-v.s).trunc(d);
354         return 2;
355     }
356 
357     int getcircle(Point a,Point b,double r1,circle &c1,circle &c2){///得到过a，b两点，半径为r1的两个圆
358         circle x(a,r1),y(b,r1);
359         int t=x.pointcrosscircle(y,c1.p,c2.p);
360         if(!t) return 0;
361         c1.r=c2.r=r;
362         return t;
363     }
364 
365     int getcircle(Line u,Point q,double r1,circle &c1,circle &c2){///得到与直线u相切，过点q，半径为r1的圆
366         double dis = u.dispointtoline(q);
367         if(sgn(dis-r1*2)>0) return 0;
368         if(sgn(dis)==0) {
369             c1.p=q+((u.e-u.s).rotleft().trunc(r1));
370             c2.p=q+((u.e-u.s).rotright().trunc(r1));
371             c1.r=c2.r=r1;
372             return 2;
373         }
374         Line u1=Line((u.s+(u.e-u.s).rotleft().trunc(r1)),(u.e+(u.e-u.s).rotleft().trunc(r1)));
375         Line u2=Line((u.s+(u.e-u.s).rotright().trunc(r1)),(u.e+(u.e-u.s).rotright().trunc(r1)));
376         circle cc=circle(q,r1);
377         Point p1,p2;
378         if(!cc.pointcrossline(u1,p1,p2)) cc.pointcrossline(u2,p1,p2);
379         c1=circle(p1,r1);
380         if(p1==p2){
381             c2=c1;
382             return 1;
383         }
384         c2=circle(p2,r1);
385         return 2;
386     }
387 
388     int getcircle(Line u,Line v,double r1,circle &c1,circle &c2,circle &c3,circle &c4){///同时与直线u，v相切，半径为r1的圆
389         if(u.parallel(v)) return 0;///两直线平行
390         Line u1=Line(u.s+(u.e-u.s).rotleft().trunc(r1),u.e+(u.e-u.s).rotleft().trunc(r1));
391         Line u2=Line(u.s+(u.e-u.s).rotright().trunc(r1),u.e+(u.e-u.s).rotright().trunc(r1));
392         Line v1=Line(v.s+(v.e-v.s).rotleft().trunc(r1),v.e+(v.e-v.s).rotleft().trunc(r1));
393         Line v2=Line(v.s+(v.e-v.s).rotright().trunc(r1),v.e+(v.e-v.s).rotright().trunc(r1));
394         c1.r=c2.r=c3.r=c4.r=r1;
395         c1.p=u1.crosspoint(v1);
396         c2.p=u1.crosspoint(v2);
397         c3.p=u2.crosspoint(v1);
398         c4.p=u2.crosspoint(v2);
399         return 4;
400     }
401 
402     int getcircle(circle cx,circle cy,double r1,circle &c1,circle &c2){///同时与不相交圆 cx,cy相切，半径为r1的圆
403         circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);
404         int t=x.pointcrosscircle(y,c1.p,c2.p);
405         if(!t) return 0;
406         c1.r=c2.r=r1;
407         return t;
408     }
409 
410     int tangentline(Point q,Line &u,Line &v){///过一点作圆的切线（先判断点和圆的关系）
411         int x=relation(q);
412         if(x==2) return 0;
413         if(x==1){
414             u=Line(q,q+(q-p).rotleft());
415             v=u;
416             return 1;
417         }
418         double d=p.distance(q);
419         double l=r*r/d;
420         double h=sqrt(r*r-l*l);
421         u=Line(q,p+((q-p).trunc(l)+(q-p).rotleft().trunc(h)));
422         v=Line(q,p+((q-p).trunc(l)+(q-p).rotright().trunc(h)));
423         return 2;
424     }
425 
426     double areacircle(circle v){///求两圆相交的面积
427         int rel=relationcircle(v);
428         if(rel >= 4) return 0.0;
429         if(rel <= 2) return min(area(),v.area());
430         double d=p.distance(v.p);
431         double hf=(r+v.r+d)/2.0;
432         double ss=2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));
433         double a1=acos((r*r+d*d-v.r*v.r)/(2.0*r*d));
434         a1=a1*r*r;
435         double a2=acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));
436         a2=a2*v.r*v.r;
437         return a1+a2-ss;
438     }
439 
440     double areatriangle(Point a,Point b){///求圆和三角形pab的相交面积
441         if(sgn((p-a)^(p-b))==0) return 0.0;
442         Point q[5];
443         int len=0;
444         q[len++]=a;
445         Line l(a,b);
446         Point p1,p2;
447         if(pointcrossline(l,q[1],q[2])==2){
448             if(sgn((a-q[1])*(b-q[1]))<0) q[len++]=q[1];
449             if(sgn((a-q[2])*(b-q[2]))<0) q[len++]=q[2];
450         }
451         q[len++]=b;
452         if(len==4 && sgn((q[0]-q[1])*(q[2]-q[1]))>0) swap(q[1],q[2]);
453         double res=0;
454         for(int i=0;i<len-1;i++){
455             if(relation(q[i])==0||relation(q[i+1])==0){
456                 double arg=p.rad(q[i],q[i+1]);
457                 res+=r*r*arg/2.0;
458             }
459             else{
460                 res+=fabs((q[i]-p)^(q[i+1]-p))/2.0;
461             }
462         }
463         return res;
464     }
465 };
466 
467 struct polygon{
468     int n;
469     Point p[100010];
470     Line l[100010];
471     void input(int _n){
472         n=_n;
473         for(int i=0;i<n;i++){
474             p[i].input();
475         }
476     }
477 
478     void add(Point q){
479         p[n++]=q;
480     }
481 
482     void getline(){
483         for(int i=0;i<n;i++){
484             l[i]=Line(p[i],p[(i+1)%n]);
485         }
486     }
487 
488     struct cmp{
489         Point p;
490         cmp(const Point &p0){p=p0;}
491         bool operator()(const Point &aa,const Point &bb){
492             Point a=aa,b=bb;
493             int d=sgn((a-p)^(b-p));
494             if(d==0){
495                 return sgn(a.distance(p)-b.distance(p))<0;
496             }
497             return d>0;
498         }
499     };
500 
501     void norm(){///进行极角排序
502         Point mi=p[0];
503         for(int i=1;i<n;i++){
504             mi=min(mi,p[i]);
505         }
506         sort(p,p+n,cmp(mi));
507     }
508 
509     void getconvex(polygon &convex){///得到第一种凸包的方法，编号为0～n-1，可能要特判所有点共点或共线的特殊情况
510         sort(p,p+n);
511         convex.n=n;
512         for(int i=0;i<min(n,2);i++){
513             convex.p[i]=p[i];
514         }
515         if(convex.n==2&&(convex.p[0]==convex.p[1])) convex.n--;
516         if(n<=2) return;
517         int &top=convex.n;
518         top=1;
519         for(int i=2;i<n;i++){
520             while(top&&sgn((convex.p[top]-p[i])^(convex.p[top-1]-p[i]))<0) top--;
521             convex.p[++top]=p[i];
522         }
523         int temp=top;
524         convex.p[++top]=p[n-2];
525         for(int i=n-3;i>=0;i--){
526             while(top!=temp&&sgn((convex.p[top]-p[i])^(convex.p[top-1]-p[i]))<=0) top--;
527             convex.p[++top]=p[i];
528         }
529         if(convex.n==2&&(convex.p[0]==convex.p[1])) convex.n--;
530         convex.norm();
531     }
532 
533     void Graham(polygon &convex){///得到凸包的第二种方法
534         norm();
535         int &top=convex.n;
536         top=0;
537         if(n==1){
538             top=1;
539             convex.p[0]=p[0];
540             return;
541         }
542         if(n==2){
543             top=2;
544             convex.p[0]=p[0];
545             convex.p[1]=p[1];
546             if(convex.p[0]==convex.p[1]) top--;
547             return;
548         }
549         convex.p[0]=p[0];
550         convex.p[1]=p[1];
551         top=2;
552         for(int i=2;i<n;i++){
553             while(top>1&&sgn((convex.p[top-1]-convex.p[top-2])^(p[i]-convex.p[top-2]))<=0) top--;
554             convex.p[top++]=p[i];
555         }
556         if(convex.n==2 && (convex.p[0]==convex.p[1])) convex.n--;
557     }
558 
559     bool inconvex(){///判断是不是凸的
560         bool s[3];
561         memset(s,false,sizeof(s));
562         for(int i=0;i<n;i++){
563             int j=(i+1)%n;
564             int k=(j+1)%n;
565             s[sgn((p[j]-p[i])^(p[k]-p[i]))+1]=true;
566             if(s[0]&&s[2]) return false;
567         }
568         return true;
569     }
570 
571     int relationpoint(Point q){///判断点和任意多边形的关系  3点上 2边上 1内部 0外部
572         for(int i=0;i<n;i++){
573             if(p[i]==q) return 3;
574         }
575         getline();
576         for(int i=0;i<n;i++){
577             if(l[i].pointonseg(q)) return 2;
578         }
579         int cnt=0;
580         for(int i=0;i<n;i++){
581             int j=(i+1)%n;
582             int k=sgn((q-p[j])^(p[i]-p[j]));
583             int u=sgn(p[i].y-q.y);
584             int v=sgn(p[j].y-q.y);
585             if(k>0&&u<0&&v>=0) cnt++;
586             if(k<0&&v<0&&u>=0) cnt--;
587         }
588         return cnt!=0;
589     }
590 
591     void convexcnt(Line u,polygon &po){///直线u切割凸多边形左侧  注意直线方向
592         int &top=po.n;
593         top=0;
594         for(int i=0;i<n;i++){
595             int d1=sgn((u.e-u.s)^(p[i]-u.s));
596             int d2=sgn((u.e-u.s)^(p[(i+1)%n]-u.s));
597             if(d1>=0) po.p[top++]=p[i];
598             if(d1*d2<0)po.p[top++]=u.crosspoint(Line(p[i],p[(i+1)%n]));
599         }
600     }
601 
602     double getcircumference(){///得到周长
603         double sum=0;
604         for(int i=0;i<n;i++){
605             sum+=p[i].distance(p[(i+1)%n]);
606         }
607         return sum;
608     }
609 
610     double getarea(){///得到面积
611         double sum=0;
612         for(int i=0;i<n;i++){
613             sum+=(p[i]^p[(i+1)%n]);
614         }
615         return fabs(sum)/2;
616     }
617 
618     bool getdir(){///得到方向 1表示逆时针  0表示顺时针
619         double sum=0;
620         for(int i=0;i<n;i++){
621             sum+=(p[i]^p[(i+1)%n]);
622         }
623         if(sgn(sum)>0) return 1;
624         return 0;
625     }
626 
627     Point getbarycentre(){///得到重心
628         Point ret(0,0);
629         double area=0;
630         for(int i=1;i<n-1;i++){
631             double tmp=(p[i]-p[0])^(p[i+1]-p[0]);
632             if(sgn(tmp)==0) continue;
633             area+=tmp;
634             ret.x+=(p[0].x+p[i].x+p[i+1].x)/3*tmp;
635             ret.y+=(p[0].y+p[i].y+p[i+1].y)/3*tmp;
636         }
637         if(sgn(area)) ret =ret/area;
638         return ret;
639     }
640 
641     double areacircle(circle c){///多边形和圆交的面积
642         double ans=0;
643         for(int i=0;i<n;i++){
644             int j=(i+1)%n;
645             if(sgn((p[j]-c.p)^(p[i]-c.p))>=0) ans+=c.areatriangle(p[i],p[j]);
646             else ans-=c.areatriangle(p[i],p[j]);
647         }
648         return fabs(ans);
649     }
650 
651     int relationcircle(circle c){///多边形和圆的关系 2圆完全在多边形内  1圆在多边形里面，碰到了多边形的边界  0其他
652         getline();
653         int x=2;
654         if(relationpoint(c.p)!=1) return 0;
655         for(int i=0;i<n;i++){
656             if(c.relationseg(l[i])==2) return 0;
657             if(c.relationseg(l[i])==1) x=1;
658         }
659         return x;
660     }
661 };
662 
663 double cross(Point a,Point b,Point c){///ab x ac
664     return (b-a)^(c-a);
665 }
666 
667 double dot(Point a,Point b,Point c){///ab*ac;
668     return (b-a)*(c-a);
669 }
670 
671 double minRectangleCover(polygon A){///`最小矩形面积覆盖`,` A 必须是凸包(而且是逆时针顺序)`,` 测试 UVA 10173`
672     //`要特判A.n < 3的情况`
673     if(A.n < 3)return 0.0;
674     A.p[A.n] = A.p[0];
675     double ans = -1;
676     int r = 1, p = 1, q;
677     for(int i = 0;i < A.n;i++){
678         //`卡出离边A.p[i] - A.p[i+1]最远的点`
679         while( sgn( cross(A.p[i],A.p[i+1],A.p[r+1]) - cross(A.p[i],A.p[i+1],A.p[r]) ) >= 0 )
680             r = (r+1)%A.n;
681         //`卡出A.p[i] - A.p[i+1]方向上正向n最远的点`
682         while(sgn( dot(A.p[i],A.p[i+1],A.p[p+1]) - dot(A.p[i],A.p[i+1],A.p[p]) ) >= 0 )
683             p = (p+1)%A.n;
684         if(i == 0)q = p;
685         //`卡出A.p[i] - A.p[i+1]方向上负向最远的点`
686         while(sgn(dot(A.p[i],A.p[i+1],A.p[q+1]) - dot(A.p[i],A.p[i+1],A.p[q])) <= 0)
687             q = (q+1)%A.n;
688         double d = (A.p[i] - A.p[i+1]).len2();
689         double tmp = cross(A.p[i],A.p[i+1],A.p[r]) *
690             (dot(A.p[i],A.p[i+1],A.p[p]) - dot(A.p[i],A.p[i+1],A.p[q]))/d;
691         if(ans < 0 || ans > tmp)ans = tmp;
692     }
693     return ans;
694 }
695 
696 vector<Point> convexCut(const vector<Point>&ps,Point q1,Point q2){///直线切凸多边形，多边形是逆时针的，在q1q2的左侧
697     vector<Point>qs;
698     int n=ps.size();
699     for(int i=0;i<n;i++){
700         Point p1=ps[i],p2=ps[(i+1)%n];
701         int d1=sgn((q2-q1)^(p1-q1)),d2=sgn((q2-q1)^(p2-q1));
702         if(d1>=0) qs.push_back(p1);
703         if(d1*d2<0) qs.push_back(Line(p1,p2).crosspoint(Line(q1,q2)));
704     }
705     return qs;
706 }
707 
708 struct halfplane:public Line{
709     double angle;
710     halfplane(){}
711     halfplane(Point _s,Point _e){///表示向量s->e逆时针（左侧）的半平面
712         s=_s;
713         e=_e;
714     }
715     halfplane(Line v){
716         s=v.s;
717         e=v.e;
718     }
719     void calcangle(){
720         angle=atan2(e.y-s.y,e.x-s.x);
721     }
722     bool operator<(const halfplane &b)const{
723         return angle<b.angle;
724     }
725 };
726 
727 struct halfplanes{
728     int n;
729     halfplane hp[2020];
730     Point p[2020];
731     int que[2020];
732     int st,ed;
733     void push(halfplane tmp){
734         hp[n++]=tmp;
735     }
736 
737     void unique(){///去重
738         int m=1;
739         for(int i=1;i<n;i++){
740             if(sgn(hp[i].angle-hp[i-1].angle)!=0) hp[m++]=hp[i];
741             else if(sgn((hp[m-1].e-hp[m-1].s)^(hp[i].s-hp[m-1].s))>0) hp[m-1]=hp[i];
742         }
743         n=m;
744     }
745     bool halfplaneinsert(){
746         for(int i=0;i<n;i++) hp[i].calcangle();
747         sort(hp,hp+n);
748         unique();
749         que[st=0]=0;
750         que[ed=1]=1;
751         p[1]=hp[0].crosspoint(hp[1]);
752         for(int i=2;i<n;i++){
753             while(st<ed&&sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<0) ed--;
754             while(st<ed&&sgn((hp[i].e-hp[i].s)^(p[st+1]-hp[i].s))<0) st++;
755             que[++ed]=i;
756             if(hp[i].parallel(hp[que[ed-1]])) return false;
757             p[ed]=hp[i].crosspoint(hp[que[ed-1]]);
758         }
759         while(st<ed&&sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<0) ed--;
760         while(st<ed&&sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+1]-hp[que[ed]].s))<0) st++;
761         if(st+1>=ed) return false;
762         return true;
763     }
764 
765     void getconvex(polygon &con){///得到最后半平面交得到的凸多边形，要先调用halfplaneinsert()且返回true
766         p[st]=hp[que[st]].crosspoint(hp[que[ed]]);
767         con.n=ed-st+1;
768         for(int j=st,i=0;j<=ed;i++,j++){
769             con.p[i]=p[j];
770         }
771     }
772 };
773 
774 struct circles{
775     circle c[1010];
776     double ans[1010];///ans[i]表示被覆盖了i次的面积
777     double pre[1010];
778     int n;
779     circles(){}
780     void add(circle cc){
781         c[n++]=cc;
782     }
783 
784     bool inner(circle x,circle y){///x包含在y中
785         if(x.relationcircle(y)!=1) return 0;
786         return sgn(x.r-y.r)<=0?1:0;
787     }
788 
789     void init_or(){///圆的面积并去掉内含的圆
790         bool mark[1010]={0};
791         int i,j,k=0;
792         for(i=0;i<n;i++){
793             for(j=0;j<n;j++){
794                 if(i!=j&&!mark[j]){
795                     if(c[i]==c[j]||inner(c[i],c[j])) break;
796                 }
797             }
798             if(j<n) mark[i]=1;
799         }
800         for(i=0;i<n;i++){
801             if(!mark[i]) c[k++]=c[i];
802         }
803         n=k;
804     }
805 
806     void init_add(){///圆的面积交去掉内含的圆
807         int i,j,k;
808         bool mark[1010]={0};
809         for(i=0;i<n;i++){
810             for(int j=0;j<n;j++){
811                 if(i!=j&&!mark[j]){
812                     if((c[i]==c[j])||inner(c[j],c[i])) break;
813                 }
814             }
815             if(j<n) mark[i]=1;
816         }
817         for(i=0;i<n;i++){
818             if(!mark[i]){
819                 c[k++]=c[i];
820             }
821         }
822         n=k;
823     }
824 
825     double areaarc(double th,double r){///半径为r的圆，弧度为th，对应的弓形的面积
826         return 0.5*r*r*(th-sin(th));
827     }
828 
829 
830     void getarea(){
831         memset(ans,0,sizeof(ans));
832         vector<pair<double,int> >v;
833         for(int i=0;i<n;i++){
834             v.clear();
835             v.push_back(make_pair(-PI,1));
836             v.push_back(make_pair(PI,-1));
837             for(int j=0;j<n;j++){
838                 if(i!=j){
839                     Point q=(c[j].p-c[i].p);
840                     double ab=q.len(),ac=c[i].r,bc=c[j].r;
841                     if(sgn(ab+ac-bc)<=0){
842                         v.push_back(make_pair(-PI,1));
843                         v.push_back(make_pair(PI,-1));
844                         continue;
845                     }
846                     if(sgn(ab+bc-ac)<=0)continue;
847                     if(sgn(ab-ac-bc)>0) continue;
848                     double th=atan2(q.y,q.x),fai=acos((ac*ac+ab*ab-bc*bc)/(2.0*ac*ab));
849                     double a0=th-fai;
850                     if(sgn(a0+PI)<0) a0+=2*PI;
851                     double a1=th+fai;
852                     if(sgn(a1-PI)>0) a1-=2*PI;
853                     if(sgn(a0-a1)>0){
854                         v.push_back(make_pair(a0,1));
855                         v.push_back(make_pair(PI,-1));
856                         v.push_back(make_pair(-PI,1));
857                         v.push_back(make_pair(a1,-1));
858                     }
859                     else{
860                         v.push_back(make_pair(a0,1));
861                         v.push_back(make_pair(a1,-1));
862                     }
863                 }
864                 sort(v.begin(),v.end());
865                 int cur=0;
866                 for(int j=0;j<v.size();j++){
867                     if(cur&&sgn(v[j].first-pre[cur])){
868                         ans[cur]+=areaarc(v[j].first-pre[cur],c[i].r);
869                         ans[cur]+=0.5*(Point(c[i].p.x+c[i].r*cos(pre[cur]),c[i].p.y+c[i].r*sin(pre[cur]))^Point(c[i].p.x+c[i].r*cos(v[j].first),c[i].p.y+c[i].r*sin(v[j].first)));
870                     }
871                     cur+=v[j].second;
872                     pre[cur]=v[j].first;
873                 }
874             }
875         }
876         for(int i=1;i<n;i++){
877             ans[i]-=ans[i+1];
878         }
879     }
880 };
881 
882 
883 bool Check(Line a,Line b){
884     if(sgn((a.s-a.e)^(b.s-a.e))*sgn((a.s-a.e)^(b.e-a.e))>0) return false;
885     if(sgn((b.s-b.e)^(a.s-b.e))*sgn((b.s-b.e)^(a.e-b.e))>0) return false;
886     if(sgn(max(a.s.x,a.e.x)-min(b.s.x,b.e.x))>=0&&sgn(max(b.s.x,b.e.x)-min(a.s.x,a.e.x))>=0
887      &&sgn(max(a.s.y,a.e.y)-min(b.s.y,b.e.y))>=0&&sgn(max(b.s.y,b.e.y)-min(a.s.y,a.e.y))>=0)
888         return true;
889     else return false;
890 }
891 
892 halfplanes hps;
893 int n;
894 Point p[105];
895 polygon poly;
896 
897 int main(){
898     int t;
899     scanf("%d",&t);
900     while(t--){
901         int n;
902         scanf("%d",&n);
903         for(int i=0;i<n;i++){
904             p[i].input();
905         }
906         hps.n=0;
907         halfplane hp;
908         for(int i=0;i<n;i++){
909             hp.e=p[i];
910             hp.s=p[(i+1)%n];
911             hps.push(hp);
912         }
913         if(hps.halfplaneinsert()){
914             hps.getconvex(poly);
915             printf("%.2f\n",poly.getarea());
916         }
917         else{
918             printf("0.00\n");
919         }
920 
921     }
922     return 0;
923 }
View Code

 

posted on 2019-05-05 18:27  Fighting_sh  阅读(564)  评论(0)  编辑  收藏  举报

导航