计算几何之半平面交算法模板及应用
自己在学习半平面交的时候看到网上有几张图片很不错就 down 下来在本文中使用了……如果博主大大不乐意请联系我删除。下面给出两个我觉得不错的博客：
1.半平面交讲解1
2.半平面交讲解2
定义：
半平面：顾名思义，就是平面的一半。一条直线会把平面分成两部分，就是两个半平面。对于半平面，我们可以用直线方程式如：ax + by >= c 表示，更常用的是用向量表示，默认向量的左侧为我们所需要的半平面。
半平面交：顾名思义，就是多个半平面求交集。其结果可能是一个凸多边形、无穷平面、直线、线段、点等。
多边形的核：如果多边形中存在一个区域使得在区域中可以看到多边形中任意位置（反之亦然），则这个区域就是多边形的核。可以用半平面交来求解。
极点：听说过极坐标吧？极点就是点 (x,y) 与原点的连线与 x 轴的夹角，其范围为 [0,360].
应用：
1.判断多边形的核是否存在。
2.求多边形中可以放入的最大圆半径。
代码：
下面给出两种不同风格的代码，一种是网上找的，便于理解，一种是从kuangbin模板上摘下来的，适用性更广，可以当作板子用。其中对于给出点的顺时针和逆时针顺序不同，邝斌模板只需要加个 reverse 函数将点的顺序颠倒，而第一个模板需要更改 cmp函数和 Judge函数中的 < 或 >。
int n,tol,bot,top; //bot和top指向双端队列的底部和顶部
int dq[MAXN],order[MAXN]; //双端队列和存储排序后点的下标
struct point
{ //点
double x,y;
} p[MAXN];
struct vct
{ //向量（直线）
point start,end;
double angle; //极角
} v[MAXN],tmp[MAXN];
int dbcmp(double k)
{ //判断为0，还是为正负
if(fabs(k)<eps) return 0;
return k>0? 1:-1;
}
double Multi(point p0,point p1,point p2)
{ //向量 p0p1 和 p0p2 相乘
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
int cmp(int x,int y)
{ //排序的比较函数
int d=dbcmp(v[x].angle-v[y].angle);
if(!d) return dbcmp(Multi(v[x].start,v[y].start,v[y].end))>0;
//大于0取向量左半部分为半平面，小于0，取右半部分
return d<0;
}
void addline(double x1,double y1,double x2,double y2)
{ //增加一个向量（直线）
v[tol].start.x=x1;
v[tol].start.y=y1;
v[tol].end.x=x2;
v[tol].end.y=y2;
//atan2是求原点到(x,y)的直线与x轴的夹角，比atan稳定
v[tol].angle=atan2(y2-y1,x2-x1);
//order[tol]=tol;
tol++;
}
void GetIntersect(vct v1,vct v2,point &p)
{ //获取交点并返回给 p
double dot1,dot2;
dot1=Multi(v1.start,v2.end,v1.end);
dot2=Multi(v1.start,v2.start,v1.end);
p.x=(dot1*v2.start.x-dot2*v2.end.x)/(dot1-dot2);
p.y=(dot1*v2.start.y-dot2*v2.end.y)/(dot1-dot2);
}
int Judge(vct v0,vct v1,vct v2)
{ //判断 v1 和 v2 的交点 pt 在半平面 v0内
point pt;
GetIntersect(v1,v2,pt);
return dbcmp(Multi(pt,v0.start,v0.end))<0;
//大于0 pt在 v0 的左面，小于0 pt 在 v0 的右面
//如果 pt 在半平面 v0 内则返回真
}
int HalfPlaneIntersection()
{ //求解半平面交
int i,j;
for(i=0;i<tol;i++) order[i]=i;
sort(order,order+n,cmp); //极角排序
for(i=1,j=0;i<tol;i++) //排除极角相同的
if(dbcmp(tmp[order[i]].angle-tmp[order[j]].angle)>0)
order[++j]=order[i];
n=j+1; //个数
dq[0]=order[0]; //开始时入队两条直线
dq[1]=order[1];
bot=0;
top=1;
for(i=2;i<tol;i++)
{ //如果双端队列底端或顶端两个半平面的交点在当前半平面之外，则删除
while(bot<top&&Judge(tmp[order[i]],tmp[dq[top-1]],tmp[dq[top]]))
top--;
while(bot<top&&Judge(tmp[order[i]],tmp[dq[bot+1]],tmp[dq[bot]]))
bot++;
dq[++top]=order[i]; //将当前半平面加入队列的顶端
}
while(bot<top&&Judge(tmp[dq[bot]],tmp[dq[top-1]],tmp[dq[top]])) top--;
while(bot<top&&Judge(tmp[dq[top]],tmp[dq[bot+1]],tmp[dq[bot]])) bot++;
if(top-bot<=1) return 0;
else return 1;
}
double GetDis(point a,point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
void ChangePolygon(double h)
{ //将所有边往内缩进 h距离
int i;
double dx,dy,len;
for(i=0;i<tol;i++)
{
len=GetDis(v[i].start,v[i].end);
dx=(v[i].start.y-v[i].end.y)/len*h;
dy=(v[i].end.x-v[i].start.x)/len*h;
tmp[i].start.x=v[i].start.x+dx;
tmp[i].start.y=v[i].start.y+dy;
tmp[i].end.x=v[i].end.x+dx;
tmp[i].end.y=v[i].end.y+dy;
tmp[i].angle=v[i].angle;
}
}
double BSearch()
{ //二分搜索
double l=0,r=20000,mid;
while(l+eps<r)
{
mid=(l+r)/2;
ChangePolygon(mid); //将所有边往内缩进 mid距离
if(HalfPlaneIntersection()) l=mid; //如果存在半平面交
else r=mid;
}
return l;
}
int sgn(double x)
{ //符号函数
if(fabs(x) < eps) return 0;
if(x < 0) return -1;
else return 1;
}
struct Point
{ //点
double x,y;
Point(){}
Point(double _x,double _y)
{
x = _x; y = _y;
}
Point operator -(const Point &b)const
{
return Point(x - b.x, y - b.y);
}
double operator ^(const Point &b)const
{ //叉积
return x*b.y - y*b.x;
}
double operator *(const Point &b)const
{ //点积
return x*b.x + y*b.y;
}
};
struct Line
{ //向量
Point s,e; //两点
double k; //斜率
Line(){}
Line(Point _s,Point _e)
{ //构造
s = _s; e = _e;
k = atan2(e.y - s.y,e.x - s.x);
}
Point operator &(const Line &b)const
{ //求两直线交点
Point res = s;
double t = ((s - b.s)^(b.s - b.e))/((s - e)^(b.s - b.e));
res.x += (e.x - s.x)*t;
res.y += (e.y - s.y)*t;
return res;
}
};
Line Q[MAXN];
Point p[MAXN]; //记录最初给的点集
Line line[MAXN]; //由最初的点集生成直线的集合
Point pp[MAXN]; //记录半平面交的结果的点集
//半平面交，直线的左边代表有效区域
bool HPIcmp(Line a,Line b)
{ //直线排序函数
if(fabs(a.k - b.k) > eps)return a.k < b.k; //斜率排序
//斜率相同我也不知道怎么办
return ((a.s - b.s)^(b.e - b.s)) < 0;
}
void HPI(Line line[], int n, Point res[], int &resn)
{ //line是半平面交的直线的集合 n是直线的条数 res是结果
//的点集 resn是点集里面点的个数
int tot = n;
sort(line,line+n,HPIcmp);
tot = 1;
for(int i = 1;i < n;i++)
if(fabs(line[i].k - line[i-1].k) > eps) //去掉斜率重复的
line[tot++] = line[i];
int head = 0, tail = 1;
Q[0] = line[0];
Q[1] = line[1];
resn = 0;
for(int i = 2; i < tot; i++)
{
if(fabs((Q[tail].e-Q[tail].s)^(Q[tail-1].e-Q[tail-1].s)) < eps ||
fabs((Q[head].e-Q[head].s)^(Q[head+1].e-Q[head+1].s)) < eps)
return;
while(head < tail && (((Q[tail]&Q[tail-1]) -
line[i].s)^(line[i].e-line[i].s)) > eps)
tail--;
while(head < tail && (((Q[head]&Q[head+1]) -
line[i].s)^(line[i].e-line[i].s)) > eps)
head++;
Q[++tail] = line[i];
}
while(head < tail && (((Q[tail]&Q[tail-1]) -
Q[head].s)^(Q[head].e-Q[head].s)) > eps)
tail--;
while(head < tail && (((Q[head]&Q[head-1]) -
Q[tail].s)^(Q[tail].e-Q[tail].e)) > eps)
head++;
if(tail <= head + 1) return;
for(int i = head; i < tail; i++)
res[resn++] = Q[i]&Q[i+1];
if(head < tail - 1)
res[resn++] = Q[head]&Q[tail];
}
double dist(Point a,Point b)
{ //两点间距离
return sqrt((a-b)*(a-b));
}
void change(Point a,Point b,Point &c,Point &d,double p)
{ //将线段ab往左移动距离p,修改得到线段cd
double len=dist(a,b);
/*三角形相似推出下面公式*/
double dx=(a.y-b.y)*p/len;
double dy=(b.x-a.x)*p/len;
c.x=a.x+dx; c.y=a.y+dy;
d.x=b.x+dx; d.y=b.y+dy;
}
double BSearch()
{ //二分搜索
double l=0,r=100000;
double ans=0;
while(r-l>=eps)
{
double mid=(l+r)/2;
for(int i=0;i < n;i++)
{
Point t1,t2;
change(p[i],p[(i+1)%n],t1,t2,mid);
line[i]=Line(t1,t2);
}
int resn;
HPI(line,n,pp,resn);
//等于0说明移多了
if(resn==0) r=mid-eps;
else l=mid+eps;
}
return l;
}半平面交算法：
最基本的是时间复杂度为 O(n) 的算法，常用的是时间复杂度为 O(nlogn) 的排序增量算法。
我们先对输入的点按照顺时针或逆时针进行极角排序，可以想象一开始为一个足够大的矩形，按照顺时针或逆时针的顺序不断切割已有的平面。求 n 个半平面的交就是用第 n 条表示当前半平面的直线去切割现有的面。每次切割都会产生一个更小的面，当所有直线都切割完毕后就是我们所需要的结果了。
排序时，我们将向量平移至以原点为起点，将坐标轴以 x 轴为界分为上下两部（x 轴属于下半部分），当两个向量终点的 y 都在 x 轴上时，按 x 从小到大排序。当两个向量重点同在上部/下部时，按叉积排序（平行按左右排序），当一上一下时，下部的排在前。
在切割的时候我们可以用一个双端队列来保存属于半平面交内的点。如上图所示，当用一条新的直线切割已有平面的时候，可能不在半平面交内的点只可能是两端的点，所以每次就不用全部点判断一次，只需要判断双端队列两端的点就可以了。
多边形的核求解方法：
前面说过判断多边形是否存在核可以用半平面交来解决。如上图所示，左上角是要求的多边形。我们可以按照顺时针的顺序用多边形的边切割一个足够大的矩形，最终剩下的就是多边形的核。
求多边形中能够放入的最大圆的半径：
可以放入的最大圆的圆心（不是圆本身）一定在多边形的核内，别问我为什么，鬼才知道……我们可以假设，最大圆半径为 r 。如果我们可以将多边形的每条边向内推进 r 长度的距离，这时候可以想见多边形里面放不开这个圆了，并且圆心所在的多边形的核也一定是变成了一个点。否则就多边形的边就还可以向内推进，这个圆就不是最大的了。
当我们不知道最大圆半径的时候，我们可以用二分的方法求半径的值，如果对于当前假设的圆半径 r ，每条边向内推进 r 距离后，发现多边形的核不存在了，则说明当前圆半径太大了，反之，就是圆半径小了，用二分找到最大的圆半径即可。
例题：
1.POJ 3335
2.POJ 3130
题意：其问题都是判断多边形的核是否存在，具体思路当然是用半平面交来解决了。要注意的是两个题给出的点的数序不同，一个是顺时针一个是逆时针，所以代码中有些地方也会有所差异。
3.POJ 3525
题意：求多边形内到边上距离最远的距离。可以转化为求多边形内可以放入的最大的圆的半径。
题解：POJ 3525 题解。
4.POJ 3384
题意：将两个半径为 r 的圆放入一个多边形中，两圆可以重叠，问两个圆占据的最大面积时的圆心坐标是多少，可以转化为上一题的方法解决。
题解：POJ 3384 题解。
推荐阅读