半平面交

半平面交可以解决一些线性规划问题和纯几何问题，比较常用于求多边形核等。
半平面交有两种写法，一种是比较直观的$O(n^2)$写法，一种是$O(n \log n)$写法。
$O(n^2)$的写法的思路是：设当前已经得到了一个半平面交的有序点集p（若半平面交封闭则必然是凸包），对于当前正在处理的直线，$O(n)$查找它能否通过切割将凸包的面积缩小，同时算出它与凸包边界的交点，从而更新凸包。
$O(n \log n)$的思想于此类似，但是使用双端队列优化，主要思路是：

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底部的半平面
step5:计算出deque顶端和底部的交点即可。

首先第一步将所有直线按照斜率排序，同时每条直线都确定一个方向使得半平面一定在它的“左侧”，这样对于斜率相同的直线一定只有最“靠左“的直线会起作用。
然后先加入前两个半平面，接着依次加入半平面。最后删除多余的不符合条件的半平面。最后要记得计算顶端和底部的交点使凸包封闭。

对于两点式直线求交的问题，可以直接列式解方程。
$y_1=k_1x+b_1$，$k_1=\frac{l_1.y_2-l_1.y_1}{l_1.x_2-l_1.x_1}$（$l_1.y_2$表示第一条直线的第二个点的纵坐标）， $y_2$同理
将上式代入$k_1x+b_1=k_2x+b_2$化简即可。

最后一份代码是点斜式求直线方程的。

用两点式时对于平行与x轴或y轴的直线也一样处理，解直线不等式时，对于直线$ax+by+c \leqslant 0$,若$a=0$则加入$(0,-c/b)$和$(1,-c/b)$，若$b=0$则加入$(-c/a,0)$和$(-c/a,1)$，若均不为0则加入$(0,-c/b)$和$(1,-a/b-c/b)$。

$O(n^2)$的写法：

POJ3335求多边形核：

 1 #include<cstdio>
 2 #include<cmath>
 3 #include<algorithm>
 4 #define rep(i,l,r) for (i=l; i<=r; i++)
 5 using namespace std;
 6 
 7 const int N=210,inf=0x3f3f3f3f;
 8 const double eps=1e-8;
 9 inline double sgn(double x){ return (fabs(x)<eps) ? 0 : ( x>0 ? 1 : -1); }
10 
11 int i,T,ntot,n,tot;
12 struct P{ double x,y; }p[N],a[N],s[N];
13 struct S{ P s,t; };
14 
15 double cross(P a,P b,P c) { return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x); }
16 bool outside(S seg,P p){ return cross(seg.s,seg.t,p)>eps; }
17 bool inside(S seg,P p){ return cross(seg.s,seg.t,p)<-eps; }
18 
19 void work(P p1,P p2,P p3,P p4){
20     double t1=cross(p1,p2,p3),t2=cross(p1,p2,p4);
21     if (sgn(t1)*sgn(t2)<0){
22         double x=(fabs(t2)*p3.x+fabs(t1)*p4.x)/(fabs(t1)+fabs(t2));
23         double y=(fabs(t2)*p3.y+fabs(t1)*p4.y)/(fabs(t1)+fabs(t2));
24         a[++tot]=(P){x,y};
25     }
26 }
27 
28 void cut(S seg){
29     int i; tot=0;
30     rep(i,1,ntot){
31         if (!outside(seg,p[i])) a[++tot]=p[i];
32         else work(seg.s,seg.t,p[i-1],p[i]),work(seg.s,seg.t,p[i],p[i+1]);
33     }
34     ntot=tot; swap(a,p); p[0]=p[tot]; p[tot+1]=p[1];
35 }
36 
37 void solve(){
38     rep(i,1,n) scanf("%lf%lf",&s[i].x,&s[i].y); s[n+1]=s[1];
39     p[1]=(P){-inf,-inf}; p[2]=(P){-inf,inf}; p[3]=(P){inf,inf}; p[4]=(P){inf,-inf};
40     ntot=4; p[0]=p[4]; p[5]=p[1];
41     rep(i,1,n) cut((S){s[i],s[i+1]});
42     if (tot==0) puts("NO"); else puts("YES");
43 }
44 
45 int main(){
46     for (scanf("%d",&T); T--; ) scanf("%d",&n),solve();
47     return 0;
48 }

POJ2540(求直线不等式的解集)

 1 #include<cstdio>
 2 #include<cmath>
 3 #include<algorithm>
 4 #define rep(i,l,r) for (i=l; i<=r; i++)
 5 using namespace std;
 6 
 7 const double eps=1e-8;
 8 char str[10];
 9 int np,nq;
10 struct P { double x,y; }pre,cur,p[110],q[110];
11 double det(P a,P b,P c){ return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y); }
12 
13 void get(P p1,P p2,double &a,double &b,double &c){
14     if (fabs(p1.x-p2.x)<eps) a=0,b=1;
15     else if (fabs(p1.y-p2.y)<eps) a=1,b=0;
16         else b=p2.y-p1.y,a=p2.x-p1.x;
17     c=-a*(p2.x+p1.x)/2-b*(p2.y+p1.y)/2;
18 }
19 
20 void work(P p1,P p2,double a,double b,double c){
21     double u=a*p1.x+b*p1.y+c,v=a*p2.x+b*p2.y+c;
22     if (u*v<-eps){
23         double x=(fabs(u)*p2.x+fabs(v)*p1.x)/(fabs(u)+fabs(v));
24         double y=(fabs(u)*p2.y+fabs(v)*p1.y)/(fabs(u)+fabs(v));
25         q[++nq]=(P){x,y};
26     }
27 }
28 
29 void cut(double a,double b,double c){
30     int i; nq=0;
31     rep(i,1,np){
32         if (a*p[i].x+b*p[i].y+c>-eps) q[++nq]=p[i];
33         else work(p[i-1],p[i],a,b,c),work(p[i],p[i+1],a,b,c);
34     }
35     swap(q,p); swap(nq,np); p[np+1]=p[1]; p[0]=p[np];
36 }
37 
38 double area(int n){
39     int i; double ar=0;
40     rep(i,2,n-1) ar+=det(p[1],p[i],p[i+1]);
41     return fabs(ar)/2.0;
42 }
43 
44 int main(){
45     p[1]=(P){0,0}; p[2]=(P){0,10}; p[3]=(P){10,10}; p[4]=(P){10,0};
46     np=4; p[0]=p[4]; p[5]=p[1];
47     double ar=1.0;
48     while (scanf("%lf%lf%s",&cur.x,&cur.y,str)!=EOF){
49         double a,b,c; get(pre,cur,a,b,c);
50         if (str[0]=='C')
51             { if (a*pre.x+b*pre.y+c<-eps) a=-a,b=-b,c=-c; }
52         else
53             if (str[0]=='H')
54                 { if (a*cur.x+b*cur.y+c<-eps) a=-a,b=-b,c=-c;}
55             else ar=0;
56         if (fabs(ar)<eps) { puts("0.00"); continue; }
57         cut(a,b,c); ar=area(np);
58         printf("%.2f\n",ar); pre=cur;
59     }
60     return 0;
61 }

$O(n \log n)$写法：

BZOJ2618:

 1 #include<cmath>
 2 #include<cstdio>
 3 #include<algorithm>
 4 #define rep(i,l,r) for (int i=l; i<=r; i++)
 5 using namespace std;
 6 
 7 const int N=1010;
 8 const double eps=1e-8;
 9 struct P{
10     double x,y;
11     P(double xx=0,double yy=0):x(xx),y(yy){}
12 }p[N],a[N];
13 struct L{P a,b; double slop;}l[1005],q[1005];
14 
15 int n,k,cnt,tot;
16 double ans;
17 
18 P operator -(P a,P b){ return P(a.x-b.x,a.y-b.y); }
19 double operator *(P a,P b){ return a.x*b.y-a.y*b.x; }
20 bool operator <(L a,L b){ return (a.slop!=b.slop) ? (a.slop-b.slop)<-eps : (a.b-a.a)*(b.b-a.a)>eps; }
21 
22 P inter(L a,L b){
23     double A1=a.b.y-a.a.y,B1=a.a.x-a.b.x,C1=(a.b.x-a.a.x)*a.a.y-(a.b.y-a.a.y)*a.a.x;
24     double A2=b.b.y-b.a.y,B2=b.a.x-b.b.x,C2=(b.b.x-b.a.x)*b.a.y-(b.b.y-b.a.y)*b.a.x;
25     return P((C2*B1-C1*B2)/(A1*B2-A2*B1),(C1*A2-C2*A1)/(A1*B2-A2*B1));
26 }
27 
28 bool jud(L a,L b,L t){ P p=inter(a,b); return (t.b-t.a)*(p-t.a)<0; }
29 
30 void work(){
31     sort(l+1,l+cnt+1);
32     int L=1,R=0; tot=0;
33     rep(i,1,cnt){
34         if (l[i].slop!=l[i-1].slop) tot++;
35         l[tot]=l[i];
36     }
37     cnt=tot; tot=0;
38     q[++R]=l[1]; q[++R]=l[2];
39     rep(i,3,cnt){
40         while (L<R && jud(q[R-1],q[R],l[i])) R--;
41         while (L<R && jud(q[L+1],q[L],l[i])) L++;
42         q[++R]=l[i];
43     }
44     while (L<R && jud(q[R-1],q[R],q[L])) R--;
45     while (L<R && jud(q[L+1],q[L],q[R])) L++;
46     q[R+1]=q[L];
47     rep(i,L,R) a[++tot]=inter(q[i],q[i+1]);
48 }
49 
50 void getans(){
51     if (tot<3) return;
52     a[++tot]=a[1];
53     rep(i,1,tot) ans+=a[i]*a[i+1];
54     ans=fabs(ans)/2;
55 }
56 
57 int main(){
58     scanf("%d",&n);
59     rep(i,1,n){
60         scanf("%d",&k);
61         rep(j,1,k) scanf("%lf%lf",&p[j].x,&p[j].y);
62         p[k+1]=p[1];
63         rep(j,1,k) l[++cnt].a=p[j],l[cnt].b=p[j+1];
64     }
65     rep(i,1,cnt) l[i].slop=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x);
66     work(); getans(); printf("%.3lf\n",ans);
67     return 0;
68 }

BZOJ4445:

 1 #include<cmath>
 2 #include<cstdio>
 3 #include<algorithm>
 4 #define A double k2=(b.b.y-b.a.y)/(b.b.x-b.a.x),b2=b.a.y-k2*b.a.x
 5 #define B double k1=(a.b.y-a.a.y)/(a.b.x-a.a.x),b1=a.a.y-k1*a.a.x
 6 #define rep(i,l,r) for (int i=l; i<=r; i++)
 7 using namespace std;
 8 
 9 const int N=100010;
10 const double eps=1e-8;
11 struct P{
12     double x,y;
13     P(double xx=0,double yy=0):x(xx),y(yy){}
14 }p[N],a[N];
15 struct L{P a,b; double slop;}l[N],q[N];
16 
17 int n,k,cnt,tot;
18 double ans1,ans2;
19 
20 P operator -(P a,P b){ return P(a.x-b.x,a.y-b.y); }
21 double operator *(P a,P b){ return a.x*b.y-a.y*b.x; }
22 bool operator <(L a,L b){ return (a.slop!=b.slop) ? (a.slop-b.slop)<-eps : (a.b-a.a)*(b.b-a.a)>eps; }
23 
24 P inter(L a,L b){
25     if (a.a.x==a.b.x){ A; return (P){a.a.x,k2*a.a.x+b2}; }
26     if (b.a.x==b.b.x){ B; return (P){b.a.x,k1*b.a.x+b1}; }
27     A; B; double x=(b2-b1)/(k1-k2),y=k1*x+b1;
28     return (P){x,y};
29 }
30 
31 bool jud(L a,L b,L t){ P p=inter(a,b); return (t.b-t.a)*(p-t.a)<-eps; }
32 
33 void work(){
34     sort(l+1,l+cnt+1);
35     int L=1,R=0; tot=0;
36     rep(i,1,cnt){
37         if (abs(l[i].slop-l[i-1].slop)>eps) tot++;
38         l[tot]=l[i];
39     }
40     cnt=tot; tot=0;
41     q[++R]=l[1]; q[++R]=l[2];
42     rep(i,3,cnt){
43         while (L<R && jud(q[R-1],q[R],l[i])) R--;
44         while (L<R && jud(q[L+1],q[L],l[i])) L++;
45         q[++R]=l[i];
46     }
47     while (L<R && jud(q[R-1],q[R],q[L])) R--;
48     while (L<R && jud(q[L+1],q[L],q[R])) L++;
49     q[R+1]=q[L];
50     rep(i,L,R) a[++tot]=inter(q[i],q[i+1]);
51 }
52 
53 void get1(){
54     ans1=0;
55     if (tot<3) return;
56     a[++tot]=a[1];
57     rep(i,1,tot) ans1+=a[i]*a[i+1];
58     ans1=fabs(ans1)/2;
59 }
60 
61 void get2(){
62     ans2=0; p[n+1]=p[1];
63     rep(i,1,n+1) ans2+=p[i]*p[i+1];
64     ans2=fabs(ans2)/2;
65 }
66 
67 void ins(int i,int j){
68     double a=p[1].y-p[i].y-p[2].y+p[j].y;
69     double b=p[2].x-p[j].x-p[1].x+p[i].x;
70     double c=p[1].x*p[2].y-p[2].x*p[1].y-p[i].x*p[j].y+p[j].x*p[i].y;
71     if (abs(a)<eps && abs(b)<eps) return;
72     if (abs(a)<eps){
73         if (b>0) l[++cnt].a=P(1,-c/b),l[cnt].b=P(0,-c/b);
74             else l[++cnt].a=P(0,-c/b),l[cnt].b=P(1,-c/b);
75     }
76     if (abs(b)<eps){
77         if (a>0) l[++cnt].a=P(-c/a,0),l[cnt].b=P(-c/a,1);
78             else l[++cnt].a=P(-c/a,1),l[cnt].b=P(-c/a,0);
79     }
80     if (abs(a)>eps && abs(b)>eps){
81         if (b>0) l[++cnt].a=P(1,-a/b-c/b),l[cnt].b=P(0,-c/b);
82             else l[++cnt].a=P(0,-c/b),l[cnt].b=P(1,-a/b-c/b);
83     }
84 }
85 
86 int main(){
87     freopen("bzoj4445.in","r",stdin);
88     freopen("bzoj4445.out","w",stdout);
89     scanf("%d",&n);
90     rep(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y);
91     rep(i,2,n) ins(i,i%n+1); l[++cnt].a=p[1],l[cnt].b=p[2];
92     rep(i,1,cnt) l[i].slop=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x);
93     work(); get1(); get2(); printf("%.4lf\n",ans1/ans2);
94     return 0;
95 }

 

posted @ 2018-02-20 12:28  HocRiser  阅读(427)  评论(0)  编辑  收藏  举报