半平面交学习笔记
半平面交学习笔记
半平面
半平面：一条直线把一个平面分成的两个平面，如图，直线\(AB\)把平面分成左（上）半平面和右（下）半平面的两个平面。
那我们如何判断是在哪一个平面呢？有两种判定方法：
-
对于直线\(AB\)，我们可以通过\(A\)、\(B\)的坐标算出其解析式，不妨设为\(ax+by+c=0\)，那么当\(C\)在左平面时，将\(C\)的坐标代入直线\(AB\)的解析式，其值大于0，即\(a\times C_x+b\times C_y+c>0\)；当\(D\)在右平面时，将\(D\)的坐标也代入直线\(AB\)的解析式，其值小于0，即\(a\times D_x+b\times D_y+c<0\)。
-
（前置知识：矢量的叉乘，可见我的矢量&凸包学习笔记）对于直线\(AB\)顺次的两个端点\(A\)、\(B\)，如果\(cross(B-A,C-A)>0\)（这里的相减操作定义为\(A-B=(A_x-B_x,A_y-B_y)\)，所以\(B-A\)就是\(B\)以\(A\)为原点的坐标，\(C-A\)就是\(C\)以\(A\)为原点的坐标），那么根据右手法则，\(C\)在左（上）半平面；如果\(cross(B-A,D-A)<0\)，那么根据右手法则，\(D\)在右（下）半平面。
半平面交
半平面交，顾名思义，就是许多半平面的交集，一般我们指的是左半平面。
我们先举道例题：
做法：首先我们注意到，某个凸多边形的每条边所在直线的半平面交就是这个凸多边形：
那么我们可以把题目给出的所有凸多边形变成对应边数的直线，然后求这些直线的半平面交。
半平面交怎么求？
首先，对于每一个凸多边形，按逆时针将每一条边建成直线。
然后对于每一条直线，以起点为原点，求出它们各自的极角，并将它们按极角排序。
排序完后，如果有两条直线极角相同，那么我们肯定只留最左侧的直线，我们就要去重。
然后将第\(1\)、\(2\)条直线放入双端队列\(deque\)，\(deque\)里存的是当前半平面交的边，因为半平面交是一个凸多边形
然后每一次加入直线到\(deque\)之前：
-
计算队尾的两条直线（也就是最陡的两条）的交点，如果它在要加入的直线的右边，那么就弹出队尾。
-
计算队首的两条直线（也就是最缓的两条）的交点，如果它在要加入的直线的右边，那么就弹出队首。
然后将新直线加入队尾。
加入完所有直线后：
-
计算队尾的两条直线（也就是最陡的两条）的交点，如果它在队首的直线的右边，那么就弹出队尾。
-
计算队首的两条直线（也就是最缓的两条）的交点，如果它在队尾的直线的右边，那么就弹出队首。
然后算出每两条直线的交点，这就是半平面交——一个凸多边形的顶点。
最后直接算出其面积就好了。（算面积方法详见我的矢量&凸包学习笔记）
最后代码如下：
#include<bits/stdc++.h>
#define eps 1e-6
#define N 510
using namespace std;
struct Point
{
double x,y;
Point(){};
Point(double a,double b)
{
x=a,y=b;
}
}p[N],e[N];
struct Line
{
Point s,t;
double d;
Line(){};
Line(Point a,Point b)
{
s=a,t=b;
}
}l[N],q[N];
Point operator + (Point a,Point b)
{
return Point{a.x+b.x,a.y+b.y};
}
Point operator - (Point a,Point b)
{
return Point{a.x-b.x,a.y-b.y};
}
Point operator * (Point a,double x)
{
return Point{a.x*x,a.y*x};
}
int n,cnt,cnt1;
double ans;
double cross(Point a,Point b)
{
return a.x*b.y-b.x*a.y;
}
double dis(Point a,Point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int compare(double x,double y)
{
if(fabs(x-y)<eps)return 0;
return x-y<0?-1:1;
}
bool cmp(Line a,Line b)
{
if(!compare(a.d,b.d))return cross(a.t-a.s,b.t-a.s)>0;
return a.d<b.d;
}
bool onright(Point a,Line l)
{
return compare(cross(l.t-l.s,a-l.s),0)>0?0:1;
}
Point intersection(Line a,Line b)
{
Point u=a.s-b.s;
Point v=a.t-a.s;
Point w=b.t-b.s;
double t=cross(w,u)/cross(v,w);
return a.s+v*t;
}
double angle(Point a)
{
return atan2(a.y,a.x);
}
void half()
{
sort(l+1,l+cnt+1,cmp);
cnt1=0;
for(int i=1;i<cnt;i++)
{
if(!compare(l[i+1].d-l[i].d,0))continue;
l[++cnt1]=l[i];
}
l[++cnt1]=l[cnt];
cnt=cnt1;
int L=1,R=0;
q[++R]=l[1],q[++R]=l[2];
for(int i=3;i<=cnt;i++)
{
while(L<R&&onright(intersection(q[R],q[R-1]),l[i]))R--;
while(L<R&&onright(intersection(q[L],q[L+1]),l[i]))L++;
q[++R]=l[i];
}
while(L<R&&onright(intersection(q[R],q[R-1]),q[L]))R--;
while(L<R&&onright(intersection(q[L],q[L+1]),q[R]))L++;
q[++R]=q[L];
cnt1=0;
for(int i=L;i<R;i++)e[++cnt1]=intersection(q[i],q[i+1]);
e[++cnt1]=e[1];
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
int m;
scanf("%d",&m);
Point last,st;
for(int j=1;j<=m;j++)
{
double x,y;
scanf("%lf%lf",&x,&y);
if(j==1)
{
last=st=Point{x,y};
continue;
}
l[++cnt]=Line{last,Point{x,y}};
last=Point{x,y};
}
l[++cnt]=Line{last,st};
}
for(int i=1;i<=cnt;i++)l[i].d=angle(l[i].s-l[i].t);
half();
if(cnt1<3)
{
puts("0.000");
return 0;
}
for(int i=1;i<cnt1;i++)ans+=cross(e[i],e[i+1]);
printf("%.3lf\n",ans/2.0);
return 0;
}
练习
poj3130 How I Mathematician Wonder What You Are!/poj3335 Rotating Scoreboard/poj1474 Video Surveillance 判断半平面交是否是平面而不是一条线、一个点或没有。看最后交点数是否大于2不就好了……
代码：
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define eps 1e-6
#define N 510
using namespace std;
struct Point
{
double x,y;
Point(){};
Point(double a,double b)
{
x=a,y=b;
}
}p[N],e[N];
struct Line
{
Point s,t;
double d;
Line(){};
Line(Point a,Point b)
{
s=a,t=b;
}
}l[N],q[N];
Point operator + (Point a,Point b)
{
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b)
{
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double x)
{
return Point(a.x*x,a.y*x);
}
int n,cnt,cnt1;
double ans;
double cross(Point a,Point b)
{
return a.x*b.y-b.x*a.y;
}
double dis(Point a,Point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int compare(double x,double y)
{
if(fabs(x-y)<eps)return 0;
return x-y<0?-1:1;
}
bool cmp(Line a,Line b)
{
if(!compare(a.d,b.d))return cross(a.t-a.s,b.t-a.s)>0;
return a.d<b.d;
}
bool onright(Point a,Line l)
{
return compare(cross(l.t-l.s,a-l.s),0)>0?0:1;
}
Point intersection(Line a,Line b)
{
Point u=a.s-b.s;
Point v=a.t-a.s;
Point w=b.t-b.s;
double t=cross(w,u)/cross(v,w);
return a.s+v*t;
}
double angle(Point a)
{
return atan2(a.y,a.x);
}
void half()
{
sort(l+1,l+cnt+1,cmp);
cnt1=0;
for(int i=1;i<cnt;i++)
{
if(!compare(l[i+1].d-l[i].d,0))continue;
l[++cnt1]=l[i];
}
l[++cnt1]=l[cnt];
cnt=cnt1;
int L=1,R=0;
q[++R]=l[1],q[++R]=l[2];
for(int i=3;i<=cnt;i++)
{
while(L<R&&onright(intersection(q[R],q[R-1]),l[i]))R--;
while(L<R&&onright(intersection(q[L],q[L+1]),l[i]))L++;
q[++R]=l[i];
}
while(L<R&&onright(intersection(q[R],q[R-1]),q[L]))R--;
while(L<R&&onright(intersection(q[L],q[L+1]),q[R]))L++;
q[++R]=q[L];
cnt1=0;
for(int i=L;i<R;i++)e[++cnt1]=intersection(q[i],q[i+1]);
e[++cnt1]=e[1];
}
int main()
{
while(~scanf("%d",&n)&&n)
{
cnt=0;
Point last,st;
for(int i=1;i<=n;i++)
{
double x,y;
scanf("%lf%lf",&x,&y);
if(i==1)
{
last=st=Point(x,y);
continue;
}
l[++cnt]=Line(last,Point(x,y));
last=Point(x,y);
}
l[++cnt]=Line(last,st);
for(int i=1;i<=cnt;i++)l[i].d=angle(l[i].s-l[i].t);
half();
ans=0;
if(cnt1<3)ans=0;
else for(int i=1;i<cnt1;i++)ans+=cross(e[i],e[i+1]);
if(!ans)puts("0");
else puts("1");
}
return 0;
}
poj1279 Art Gallery 模板……代码不放了。
poj3384 Feng Shui 将凸多边形的每条边内移\(r\)，那么圆心就在这些直线的半平面交的端点上，最后暴力枚举或旋转卡壳即可。
代码：
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define eps 1e-6
#define N 110
#define INF 0x7fffffff
using namespace std;
struct Point
{
double x,y;
Point(){};
Point(double a,double b)
{
x=a,y=b;
}
}e[N],p[N],p1,p2;
Point operator + (Point a,Point b)
{
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b)
{
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double x)
{
return Point(a.x*x,a.y*x);
}
double angle(Point a)
{
return atan2(a.y,a.x);
}
struct Line
{
Point s,t;
double d;
Line(){};
Line(Point a,Point b)
{
s=a,t=b,d=angle(b-a);
}
}l[N],q[N];
int n,cnt,cnt1;
double r,ans=-INF;
double cross(Point a,Point b)
{
return a.x*b.y-b.x*a.y;
}
double dis(Point a,Point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int compare(double x,double y)
{
if(fabs(x-y)<eps)return 0;
return x-y<0?-1:1;
}
bool cmp(Line a,Line b)
{
if(!compare(a.d,b.d))return compare(cross(a.t-a.s,b.t-a.s),0)>0;
return a.d<b.d;
}
bool onright(Point a,Line l)
{
return compare(cross(l.t-l.s,a-l.s),0)>=0?0:1;
}
Point intersection(Line a,Line b)
{
Point u=a.s-b.s;
Point v=a.t-a.s;
Point w=b.t-b.s;
double t=cross(w,u)/cross(v,w);
return a.s+v*t;
}
void half()
{
sort(l+1,l+cnt+1,cmp);
cnt1=0;
for(int i=1;i<cnt;i++)
{
if(!compare(l[i+1].d-l[i].d,0))continue;
l[++cnt1]=l[i];
}
l[++cnt1]=l[cnt];
cnt=cnt1;
int L=1,R=0;
q[++R]=l[1],q[++R]=l[2];
for(int i=3;i<=cnt;i++)
{
while(L<R&&onright(intersection(q[R],q[R-1]),l[i]))R--;
while(L<R&&onright(intersection(q[L],q[L+1]),l[i]))L++;
q[++R]=l[i];
}
while(L<R&&onright(intersection(q[R],q[R-1]),q[L]))R--;
while(L<R&&onright(intersection(q[L],q[L+1]),q[R]))L++;
q[++R]=q[L];
cnt1=0;
for(int i=L;i<R;i++)e[++cnt1]=intersection(q[i],q[i+1]);
}
int main()
{
scanf("%d%lf",&n,&r);
for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
p[0]=p[n];
for(int i=n;i>=1;i--)
{
double k=r/dis(p[i-1],p[i]);
double x1=k*(p[i-1].y-p[i].y);
double y1=k*(p[i].x-p[i-1].x);
l[++cnt]=Line(Point(p[i].x-x1,p[i].y-y1),Point(p[i-1].x-x1,p[i-1].y-y1));
}
half();
for(int i=1;i<=cnt1;i++)
{
for(int j=i;j<=cnt1;j++)
{
double d=dis(e[i],e[j]);
if(d>ans)
{
ans=d;
p1=e[i];
p2=e[j];
}
}
}
printf("%.4lf %.4lf %.4lf %.4lf\n",p1.x,p1.y,p2.x,p2.y);
return 0;
}
