再来一道测半平面交模板题 Poj1279 Art Gallery

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cmath>
  4 #include <algorithm>
  5 
  6 
  7 using namespace std;
  8 const double eps = 1e-8;
  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 = 2000;
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;
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;
474         int bot = 0, top = 1;
475         dequeue[0] = _Off[0];
476         dequeue[1] = _Off[1];
477         for(i = 2; i < ln; i ++)
478         {
479             while(bot < top &&  Polygon::judege(dequeue[top],dequeue[top-1],_Off[i]) > 0)
480                 top --;
481             while(bot < top &&  Polygon::judege(dequeue[bot],dequeue[bot+1],_Off[i]) > 0)
482                 bot ++;
483             dequeue[++ top] = _Off[i];
484         }
485 
486         while(bot < top && Polygon::judege(dequeue[top],dequeue[top-1],dequeue[bot]) > 0)
487             top --;
488         while(bot < top && Polygon::judege(dequeue[bot],dequeue[bot+1],dequeue[top]) > 0)
489             bot ++;
490         //计算交点(注意不同直线形成的交点可能重合)
491         if(top <= bot + 1)
492             return (*this);
493         for(i = bot; i < top; i ++)
494             Line::crossLPt(dequeue[i],dequeue[i + 1],pt[n++]);
495         if(bot < top + 1)
496             Line::crossLPt(dequeue[bot],dequeue[top],pt[n++]);
497         return (*this);
498     }
499 };
500 
501 int n,t;
502 Point pt[2000];
503 Line ln[2000];
504 Polygon ans;
505 int main(void)
506 {
507     scanf("%d",&t);
508     while(t--)
509     {
510         scanf("%d",&n);
511         for(int i=0;i<n;i++)
512             scanf("%lf%lf",&pt[i].x,&pt[i].y);
513         pt[n++]=pt[0];
514         for(int i=n-1;i;i--)
515             ln[i-1]=Line(pt[i],pt[i-1]);
516         //for(int i=0;i<n-1;i++)
517         //    printf("%.2f %.2f %.2f %.2f\n",ln[i].s.x,ln[i].s.y,ln[i].e.x,ln[i].e.y);
518         ans.halfPanelCross(ln,n-1);
519         double area=0;
520         if(ans.n==0)
521             area=0;
522         else
523             area=ans.area();
524         printf("%.2f\n",area);
525     }
526     return 0;
527 }