计算几何学习之半平面交
程序员文章站
2022-04-01 13:32:33
...
首先解决问题：什么是半平面？ 顾名思义，半平面就是指平面的一半，我们知道，一条直线可以将平面分为两个部分，那么这两个部分就叫做两个半平面。
然后，半平面怎么表示呢？ 二维坐标系下，直线可以表示为ax + by + c = 0，那么两个半平面则可以表示为ax + by + c >= 0 和ax + by + c < 0，这就是半平面的表示方法。
还有，半平面的交是神马玩意？ 其实就是一个方程组，让你画出满足若干个式子的坐标系上的区域（类似于线性规划的可行域），方程组就是由类似于上面的这些不等式组成的。
另外，半平面交可以干什么？ 半平面交虽然说是半平面的问题，但它其实就是关于直线的问题。一个一个的半平面其实就是一个一个有方向的直线而已。
半平面交的一个重要应用就是求多边形的核 。 多边形的核又是神马玩意？ 它是平面简单多边形的核是该多边形内部的一个点集，该点集中任意一点与多边形边界上一点的连线都处于这个多边形内部。就是一个在一个房子里面放一个摄像 头，能将所有的地方监视到的放摄像头的地点的集合即为多边形的核。经常会遇到让你判定一个多边形是否有核的问题。
比如说：
这三个题比较简单，裸的半平面交，下面是核心代码（暨半平面交的模板）：
/*半平面相交（直线切割多边形）（点标号从1开始）*/
Point points[MAXN],p[MAXN],q[MAXN];
int n;
double r;
int cCnt,curCnt;
inline void getline(Point x,Point y,double &a,double &b,double &c){
a = y.y - x.y;
b = x.x - y.x;
c = y.x * x.y - x.x * y.y;
}
inline void initial(){
for(int i = 1; i <= n; ++i)p[i] = points[i];
p[n+1] = p[1];
p[0] = p[n];
cCnt = n;
}
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
inline void cut(double a,double b ,double c){
curCnt = 0;
for(int i = 1; i <= cCnt; ++i){
if(a*p[i].x + b*p[i].y + c >= EPS)q[++curCnt] = p[i];
else {
if(a*p[i-1].x + b*p[i-1].y + c > EPS){
q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
}
if(a*p[i+1].x + b*p[i+1].y + c > EPS){
q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
}
}
}
for(int i = 1; i <= curCnt; ++i)p[i] = q[i];
p[curCnt+1] = q[1];p[0] = p[curCnt];
cCnt = curCnt;
}
inline void solve(){
//注意：默认点是顺时针，如果题目不是顺时针，规整化方向
initial();
for(int i = 1; i <= n; ++i){
double a,b,c;
getline(points[i],points[i+1],a,b,c);
cut(a,b,c);
}
/*
如果要向内推进r，用该部分代替上个函数
for(int i = 1; i <= n; ++i){
Point ta, tb, tt;
tt.x = points[i+1].y - points[i].y;
tt.y = points[i].x - points[i+1].x;
double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);
tt.x = tt.x * k;
tt.y = tt.y * k;
ta.x = points[i].x + tt.x;
ta.y = points[i].y + tt.y;
tb.x = points[i+1].x + tt.x;
tb.y = points[i+1].y + tt.y;
double a,b,c;
getline(ta,tb,a,b,c);
cut(a,b,c);
}
*/
//多边形核的面积
double area = 0;
for(int i = 1; i <= curCnt; ++i)
area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;
area = fabs(area / 2.0);
//此时cCnt为最终切割得到的多边形的顶点数，p为存放顶点的数组
}
inline void GuiZhengHua(){
//规整化方向，逆时针变顺时针，顺时针变逆时针
for(int i = 1; i < (n+1)/2; i ++)
swap(points[i], points[n-i]);//头文件加iostream
}
inline void init(){
for(int i = 1; i <= n; ++i)points[i].input();
points[n+1] = points[1];
}
const double eps = 1e-8;
const double PI = acos(-1.0);
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;
}
};
//半平面交，直线的左边代表有效区域
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;
}
Line Q[110];
void HPI(Line line[], int n, Point res[], int &resn)//res[]是最后多边形的核，resn是核的顶点数，当resn==0时没有核
{
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];
}普通半平面交写法
poj1755
const double eps = 1e-18;
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;
}
};
//计算多边形面积
double CalcArea(Point p[],int n)
{
double res = 0;
for(int i = 0;i < n;i++)
res += (p[i]^p[(i+1)%n]);
return fabs(res/2);
}
//通过两点，确定直线方程
void Get_equation(Point p1,Point p2,double &a,double &b,double &c)
{
a = p2.y - p1.y;
b = p1.x - p2.x;
c = p2.x*p1.y - p1.x*p2.y;
}
//求交点
Point Intersection(Point p1,Point p2,double a,double b,double c)
{
double u = fabs(a*p1.x + b*p1.y + c);
double v = fabs(a*p2.x + b*p2.y + c);
Point t;
t.x = (p1.x*v + p2.x*u)/(u+v);
t.y = (p1.y*v + p2.y*u)/(u+v);
return t;
}
Point tp[110];
void Cut(double a,double b,double c,Point p[],int &cnt)
{
int tmp = 0;
for(int i = 1;i <= cnt;i++)
{
//当前点在左侧，逆时针的点
if(a*p[i].x + b*p[i].y + c < eps)tp[++tmp] = p[i];
else
{
if(a*p[i-1].x + b*p[i-1].y + c < -eps)
tp[++tmp] = Intersection(p[i-1],p[i],a,b,c);
if(a*p[i+1].x + b*p[i+1].y + c < -eps)
tp[++tmp] = Intersection(p[i],p[i+1],a,b,c);
}
}
for(int i = 1;i <= tmp;i++)
p[i] = tp[i];
p[0] = p[tmp];
p[tmp+1] = p[1];
cnt = tmp;
}
double V[110],U[110],W[110];
int n;
const double INF = 100000000000.0;
Point p[110];
bool solve(int id)
{
p[1] = Point(0,0);
p[2] = Point(INF,0);
p[3] = Point(INF,INF);
p[4] = Point(0,INF);
p[0] = p[4];
p[5] = p[1];
int cnt = 4;
for(int i = 0;i < n;i++)
if(i != id)
{
double a = (V[i] - V[id])/(V[i]*V[id]);
double b = (U[i] - U[id])/(U[i]*U[id]);
double c = (W[i] - W[id])/(W[i]*W[id]);
if(sgn(a) == 0 && sgn(b) == 0)
{
if(sgn(c) >= 0)return false;
else continue;
}
Cut(a,b,c,p,cnt);
}
if(sgn(CalcArea(p,cnt)) == 0)return false;
else return true;
}
int main()
{
while(scanf("%d",&n) == 1)
{
for(int i = 0;i < n;i++)
scanf("%lf%lf%lf",&V[i],&U[i],&W[i]);
for(int i = 0;i < n;i++)
{
if(solve(i))printf("Yes\n");
else printf("No\n");
}
}
return 0;
}对 凸多边形（指的是当前多边形）的每一个顶点，如果这个顶点在直线的指定的一侧（暨在该半平面上），那么就将该顶点直接加入到新的多边形中，否则，看与该点 相邻的多边形上的两个点（判断线段是否和直线相交），如果和直线相交，则把交点加入到新的多边形中。这样，就可以得到一个新的凸多边形。关于最初的多边 形，我们可以设一个四方形的区域，比如说（-1000，-1000），（-1000,1000），（1000,1000），（1000,-1000），然 后开始切割。
POJ 3525 Most Distant Point from the Sea （推荐）
求在多边形中画一个圆的最大半径
可以使用半平面交+ 二分法解。对每个距离进行判定（是否存在区域可以放置圆心，多边形的的每条边向内推进
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
using namespace std;
#define EPS 1e-10
const int MAXN = 200;
struct Point{
double x,y;
Point(){}
Point(double _x,double _y):x(_x),y(_y){}
void input(){
scanf("%lf%lf",&x,&y);
}
};
Point points[MAXN],p[MAXN],q[MAXN];
int n;
int cCnt,curCnt;
inline void getline(Point x,Point y,double &a,double &b,double &c){
a = y.y - x.y;
b = x.x - y.x;
c = y.x * x.y - x.x * y.y;
}
inline void initial(){
for(int i = 1; i <= n; ++i)p[i] = points[i];
p[n+1] = p[1];
p[0] = p[n];
cCnt = n;
}
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
inline void cut(double a,double b ,double c){
curCnt = 0;
for(int i = 1; i <= cCnt; ++i){
if(a*p[i].x + b*p[i].y + c >= 0)q[++curCnt] = p[i];
else {
if(a*p[i-1].x + b*p[i-1].y + c > 0){
q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
}
if(a*p[i+1].x + b*p[i+1].y + c > 0){
q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
}
}
}
for(int i = 1; i <= curCnt; ++i)p[i] = q[i];
p[curCnt+1] = q[1];p[0] = p[curCnt];
cCnt = curCnt;
}
inline int solve(double r){
initial();
for(int i = 1; i <= n; ++i){
Point ta, tb, tt;
tt.x = points[i+1].y - points[i].y;
tt.y = points[i].x - points[i+1].x;
double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);
tt.x = tt.x * k;
tt.y = tt.y * k;
ta.x = points[i].x + tt.x;
ta.y = points[i].y + tt.y;
tb.x = points[i+1].x + tt.x;
tb.y = points[i+1].y + tt.y;
double a,b,c;
getline(ta,tb,a,b,c);
cut(a,b,c);
}
if(cCnt <= 0)return 0;
return 1;
}
int main(){
while(scanf("%d",&n) != EOF){
if(n == 0)break;
for(int i = 1; i <= n; ++i)points[i].input();
for(int i = 1; i < (n+1)/2; i ++)
swap(points[i], points[n-i]);
points[n+1] = points[1];
double high = 10000000;
double low = 0,mid;
while(low + EPS <= high){
mid = (high + low)/2.0;
if(solve(mid))low = mid;
else high = mid;
}
printf("%lf/n",high);
}
return 0;
} 半平面交实际应用，用两个圆覆盖一个多边形，问最多能覆盖多边形的面积。
解法：用半平面交将多边形的每条边一起向“内”推进R，得到新的多边形，然后求多边形的最远两点
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
using namespace std;
#define EPS 1e-10
const int MAXN = 300;
struct Point{
double x,y;
Point(){}
Point(double _x,double _y):x(_x),y(_y){}
void input(){
scanf("%lf%lf",&x,&y);
}
};
Point points[MAXN],p[MAXN],q[MAXN];
int n;
double r;
int cCnt,curCnt;
void getline(Point x,Point y,double &a,double &b,double &c){
a = y.y - x.y;
b = x.x - y.x;
c = y.x * x.y - x.x * y.y;
}
inline void initial(){
for(int i = 1; i <= n; ++i)p[i] = points[i];
p[n+1] = p[1];
p[0] = p[n];
cCnt = n;
}
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
inline void cut(double a,double b ,double c){
curCnt = 0;
for(int i = 1; i <= cCnt; ++i){
if(a*p[i].x + b*p[i].y + c >= 0)q[++curCnt] = p[i];
else {
if(a*p[i-1].x + b*p[i-1].y + c > 0){
q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
}
if(a*p[i+1].x + b*p[i+1].y + c > 0){
q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
}
}
}
for(int i = 1; i <= curCnt; ++i)p[i] = q[i];
p[curCnt+1] = q[1];p[0] = p[curCnt];
cCnt = curCnt;
}
inline double pdis(Point a,Point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int main(){
scanf("%d%lf",&n,&r);
for(int i = 1; i <= n; ++i)points[i].input();
points[n+1] = points[1];
initial();
for(int i = 1; i <= n; ++i){
Point ta, tb, tt;
tt.x = points[i+1].y - points[i].y;
tt.y = points[i].x - points[i+1].x;
double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);
tt.x = tt.x * k;
tt.y = tt.y * k;
ta.x = points[i].x + tt.x;
ta.y = points[i].y + tt.y;
tb.x = points[i+1].x + tt.x;
tb.y = points[i+1].y + tt.y;
double a,b,c;
getline(ta,tb,a,b,c);
cut(a,b,c);
}
int ansx = 0,ansy = 0;
double res = 0;
for(int i = 1; i <= cCnt; ++i){
for(int j = i + 1; j <= cCnt; ++j){
double tmp = pdis(p[i],p[j]);
if(tmp > res){
res = tmp;
ansx = i;
ansy = j;
}
}
}
printf("%.4lf %.4lf %.4lf %.4lf/n",p[ansx].x,p[ansx].y,p[ansy].x,p[ansy].y);
return 0;
} POJ 1755 Triathlon （推荐）
半平面交判断不等式是否有解。注意特殊情况
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define EPS 1e-8
const int MAXN = 110;
struct Point{
double x,y;
Point(){}
Point(double _x,double _y):x(_x),y(_y){}
};
struct Node{
double x,y,z;
void input(){
scanf("%lf%lf%lf",&x,&y,&z);
}
};
inline void getabc(double &a,double &b,double &c,Node x1,Node x2){
a = 1.0/x2.x - 1.0/x1.x ;
b = 1.0/x2.y - 1.0/x1.y;
c = 1.0/x2.z - 1.0/x1.z;
}
Point points[MAXN],p[MAXN],q[MAXN];
int n;
int cCnt,curCnt;
int flag ;
int sig(double k){
return (k < -EPS) ? -1 : (k > EPS);
}
inline void getline(Point x,Point y,double &a,double &b,double &c){
a = y.y - x.y;
b = x.x - y.x;
c = y.x * x.y - x.x * y.y;
}
inline void initial(){
p[1] = Point(0,0);
p[2] = Point(0,1000000);
p[3] = Point(1000000,1000000);
p[4] = Point(1000000,0);
p[5] = p[1];
p[0] = p[4];
cCnt = 4;
}
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
inline void cut(double a,double b ,double c){
curCnt = 0;
if(a == 0 && b == 0 && c == 0){
flag = 0;
return ;
}
for(int i = 1; i <= cCnt; ++i){
if(sig(a*p[i].x + b*p[i].y + c) >= 0)q[++curCnt] = p[i];
//注意这里一定是大于0，而不是大于等于0，直接排除掉a = b = c = 0这种情况
else {
if(sig(a*p[i-1].x + b*p[i-1].y + c) > 0){
q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
}
if(sig(a*p[i+1].x + b*p[i+1].y + c) > 0){
q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
}
}
}
double area = 0;
for(int i = 1; i <= curCnt; ++i)
area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;
area = fabs(area / 2.0);
if(area <= EPS)flag = 0;
for(int i = 1; i <= curCnt; ++i)p[i] = q[i];
p[curCnt+1] = q[1];p[0] = p[curCnt];
cCnt = curCnt;
}
Node nodes[MAXN];
int N;
int main(){
scanf("%d",&N);
for(int i = 1; i <= N; ++i)
nodes[i].input();
bool ok;
double a,b,c;
for(int i = 1; i <= N; ++i){
ok = 1;
initial();
flag = 1;
for(int j = 1; j <= N; j++){
if(i == j)continue;
getabc(a,b,c,nodes[i],nodes[j]);
cut(a,b,c);
if(flag == 0){
ok = 0;
break;
}
}
if(ok)printf("Yes/n");
else printf("No/n");
}
return 0;
} POJ 2540 Hotter Colder
半平面交求线性规划可行区域的面积，需要求两点的中垂线。
注意直线切割多边形的时候，向新的点数组（q）里加点的时候，点的顺序必须和原来保持一致（顺时针）
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define EPS 1e-6
const int MAXN = 200;
struct Point{
double x,y;
Point(){}
Point(double _x,double _y):x(_x),y(_y){}
void input(){
scanf("%lf%lf",&x,&y);
}
};
Point p[MAXN],q[MAXN];
int cCnt ,curCnt;
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
int sig(double k){
return (k < -EPS) ? -1 : (k > EPS);
}
inline void gethalf(Point p1,Point p2,double &a,double &b,double &c,
int flag){
int fhalf;
a = p2.x - p1.x;
b = p2.y - p1.y;
c = (p1.y*p1.y - p2.y*p2.y - p2.x*p2.x + p1.x*p1.x)/2;
double tmp = a*p2.x + b*p2.y + c;
if(flag == 1){
if(sig(tmp) >= 0)fhalf = 1;
else fhalf = -1;
}
if(flag == -1){
if(sig(tmp) >= 0)fhalf = -1;
else fhalf = 1;
}
if(fhalf == -1){
a = -a;
b = -b;
c = -c;
}
}
inline void cut(double a,double b,double c){
curCnt = 0;
for(int i = 1; i <= cCnt; ++i){
if(sig(a*p[i].x + b*p[i].y + c) >= 0)q[++curCnt] = p[i];
else {
if(sig(a*p[i-1].x + b*p[i-1].y + c) > 0){
q[++curCnt] = intersect(p[i-1],p[i],a,b,c);
}
if(sig(a*p[i+1].x + b*p[i+1].y + c) > 0){
q[++curCnt] = intersect(p[i+1],p[i],a,b,c);
}
}
}
for(int i = 1; i <= curCnt; ++i)
p[i] = q[i];
p[curCnt + 1] = p[1];
p[0] = p[curCnt];
cCnt = curCnt;
}
int main(){
p[1] = Point(0,0);
p[2] = Point(10,0);
p[3] = Point(10,10);
p[4] = Point(0,10);
p[5] = p[1];
p[0] = p[4];
cCnt = 4;
double a1,b1;
char op[20];
bool mark = 1;
Point last = Point(0,0);
while(scanf("%lf%lf%s",&a1,&b1,&op) != EOF){
Point current = Point(a1,b1);
if(!mark){
puts("0.00");
continue;
}
if(!strcmp(op,"Same")){
puts("0.00");
mark = 0;
continue;
}
int fhalf = 0,flag = 0;
if(!strcmp(op,"Hotter"))flag = 1;
else if(!strcmp(op,"Colder"))flag = -1;
double a,b,c;
gethalf(last,current,a,b,c,flag);
cut(a,b,c);
double area = 0;
for(int i = 1; i <= cCnt; ++i)
area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;
area = fabs(area / 2.0);
printf("%.2lf/n",area);
last = current;
}
return 0;
} 另外还有一个数据比较多的，必须用nlogn的算法解决的半平面交的问题：http://acm.pku.edu.cn/JudgeOnline/problem?id=2451
Zzy 专为他那篇nlogn算法解决半平面交问题的论文而出的题目，至今不会自己写。