POJ 2451 Uyuw’s Concert(半平面交nlgn)

轉載請註明出處，謝謝http://blog.csdn.net/acm_cxlove/article/details/7854526      
by---cxlove

題目：給出一些向量，求出圍成的多邊形的核的面積

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

本身沒啥好說的，就是半平面交的模板題，聽說 是Zzy為了他的那篇論文專門出的題，但是資料被POJ更改過，時限很寬，n^2隨便水，而且我的nlgn比n^2耗時還多。

先引用一下nlgn演算法的介紹，在zzy論文裡也有詳細介紹

step1. 將所有半平面按極角排序，對於極角相同的，選擇性的保留一個。 O(nlogn)
step2. 使用一個雙端隊列(deque)，加入最開始2個半平面。
step3. 每次考慮一個新的半平面：
  a.while deque頂端的兩個半平面的交點在當前半平面外:刪除deque頂端的半平面
  b.while deque底部的兩個半平面的交點在當前半平面外:刪除deque底部的半平面
  c.將新半平面加入deque頂端
step4.刪除兩端多餘的半平面。
具體方法是：
a.while deque頂端的兩個半平面的交點在底部半平面外:刪除deque頂端的半平面
b.while deque底部的兩個半平面的交點在頂端半平面外:刪除deque底部的半平面
重複a,b直到不能刪除為止。

step5:計算出deque頂端和底部的交點即可。

精度我是一個個試出來的。。。最終是1e-10

#include<iostream>#include<fstream>#include<iomanip>#include<cstdio>#include<cstring>#include<algorithm>#include<cstdlib>#include<cmath>#include<set>#include<map>#include<queue>#include<stack>#include<string>#include<vector>#include<sstream>#include<cassert>#define LL long long#define eps 1e-10#define inf 10000#define zero(a) fabs(a)<eps#define N 20005using namespace std;struct Point{    double x,y;}p[N*2];struct Segment{    Point s,e;    double angle;    void get_angle(){angle=atan2(e.y-s.y,e.x-s.x);}}seg[N];int m;//叉積為正說明，p2在p0-p1的左側double xmul(Point p0,Point p1,Point p2){    return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);}Point Get_Intersect(Segment s1,Segment s2){    double u=xmul(s1.s,s1.e,s2.s),v=xmul(s1.e,s1.s,s2.e);    Point t;    t.x=(s2.s.x*v+s2.e.x*u)/(u+v);t.y=(s2.s.y*v+s2.e.y*u)/(u+v);    return t;}bool cmp(Segment s1,Segment s2){    //先按極角排序    if(s1.angle>s2.angle) return true;    //極角相等，內側的在前    else if(zero(s1.angle-s2.angle)&&xmul(s2.s,s2.e,s1.e)>-eps) return true;    return false;}void HalfPlaneIntersect(Segment seg[],int n){    sort(seg,seg+n,cmp);    int tmp=1;    for(int i=1;i<n;i++)        if(!zero(seg[i].angle-seg[tmp-1].angle))            seg[tmp++]=seg[i];    n=tmp;    Segment deq[N];    deq[0]=seg[0];deq[1]=seg[1];    int head=0,tail=1;    for(int i=2;i<n;i++){        while(head<tail&&xmul(seg[i].s,seg[i].e,Get_Intersect(deq[tail],deq[tail-1]))<-eps) tail--;        while(head<tail&&xmul(seg[i].s,seg[i].e,Get_Intersect(deq[head],deq[head+1]))<-eps) head++;        deq[++tail]=seg[i];    }    while(head<tail&&xmul(deq[head].s,deq[head].e,Get_Intersect(deq[tail],deq[tail-1]))<-eps) tail--;    while(head<tail&&xmul(deq[tail].s,deq[tail].e,Get_Intersect(deq[head],deq[head+1]))<-eps) head++;    if(head==tail) return;    m=0;    for(int i=head;i<tail;i++)        p[m++]=Get_Intersect(deq[i],deq[i+1]);    if(tail>head+1)        p[m++]=Get_Intersect(deq[head],deq[tail]);}double Get_area(Point p[],int &n){    double area=0;    for(int i=1;i<n-1;i++)        area+=xmul(p[0],p[i],p[i+1]);    return fabs(area)/2.0;}int main(){    int n;    while(scanf("%d",&n)!=EOF){        seg[0].s.x=0;seg[0].s.y=0;seg[0].e.x=10000;seg[0].e.y=0;seg[0].get_angle();        seg[1].s.x=10000;seg[1].s.y=0;seg[1].e.x=10000;seg[1].e.y=10000;seg[1].get_angle();        seg[2].s.x=10000;seg[2].s.y=10000;seg[2].e.x=0;seg[2].e.y=10000;seg[2].get_angle();        seg[3].s.x=0;seg[3].s.y=10000;seg[3].e.x=0;seg[3].e.y=0;seg[3].get_angle();        for(int i=0;i<n;i++){            scanf("%lf%lf%lf%lf",&seg[i+4].s.x,&seg[i+4].s.y,&seg[i+4].e.x,&seg[i+4].e.y);            seg[i+4].get_angle();        }        HalfPlaneIntersect(seg,n+4);        printf("%.1f\n",Get_area(p,m));    }    return 0;}