半平面交算法详解-优快云博客

法一：暴力判断每个点是否被现在半平面割。但是在强制在线的时候也是唯一的方法了。

Code：

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cctype>
#include<algorithm>
using namespace std;

int siz,cnt,n;
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); }
	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; }
	Point operator *(const double &B)const{ return Point(x*B,y*B); }
}P[55],ret[505],tmp[505];

bool ERROR=0;
Point insec(Point A,Point B,Point C,Point D)
{
	ERROR=0;
	double a=(A-C) * (A-D),b=(B-D) * (B-C);
	if(a*b<0){ ERROR=1;return Point(0,0); }
	return A+(B-A)*(a/(a+b));
}

void Cut(Point A,Point B)
{
	cnt=0;
	Point res;
	for(int i=0;i<siz;i++)
		if((A - ret[i]) * (B - ret[i]) < 0)
		{
			if((res = insec(ret[i],ret[(i-1+siz)%siz],A,B) , ERROR)==0) 
				tmp[cnt++]=res;
			if((res = insec(ret[i],ret[(i+1)%siz],A,B) , ERROR)==0) 
				tmp[cnt++]=res;
		}
		else tmp[cnt++]=ret[i];
	for(int i=0;i<cnt;i++) ret[i]=tmp[i];
	siz=cnt;
}

int main()
{
	int T;scanf("%d",&T);
	scanf("%d",&siz);
	for(int i=0;i<siz;i++) scanf("%lf%lf",&ret[i].x,&ret[i].y);
	for(T--;T--;)
	{
		scanf("%d",&n);
		for(int i=0;i<n;i++) scanf("%lf%lf",&P[i].x,&P[i].y);
		for(int i=0;i<n;i++)
			Cut(P[i],P[(i+1)%n]);
	}
	if(siz<3){ printf("0.000");return 0; }
	
	double ans=0;
	for(int i=0;i<siz;i++) 
		ans+=ret[i] * ret[(i+1)%siz] * 0.5;
	printf("%.3lf",fabs(ans));
}

法二：

nlogn半平面交。POJ 2451

AC Code:

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cctype>
#include<cmath>
#define maxn 100005
#define eps 1e-8
using namespace std;

struct Point
{
	double x,y;
	Point (double x=0,double y=0):x(x),y(y){}
	Point operator +(const Point B)const{ return Point(x+B.x,y+B.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; }
	Point operator *(const double B)const{ return Point(x*B,y*B); }
}P[maxn];

struct Line
{
	Point A,B;
	double ag;
	bool operator <(const Line &nt)const
	{   if(ag == nt.ag) return (nt.B-nt.A)*(B-nt.A)>0;
		return ag < nt.ag; }
}L[maxn];

int dcmp(double a)
{
	if(a>eps) return 1;
	if(a<-eps) return -1;
	return 0;
}
Point qp[maxn];
int q[maxn],l,r,n;

Point Ipt(Point p1,Point p2,Point p3,Point p4)
{
	Point u=p2-p1,v=p4-p3,w=p1-p3;
	return p1+(p2-p1)*((v*w)/(u*v));
}

int main()
{
	int tot;
	scanf("%d",&tot);
	for(int i=1;i<=tot;i++)
	{
		scanf("%lf%lf%lf%lf",&L[i].A.x,&L[i].A.y,&L[i].B.x,&L[i].B.y);
		L[i].ag = atan2(L[i].B.y-L[i].A.y,L[i].B.x-L[i].A.x);
	}
	L[++tot].A = Point(0,0) , L[tot].B = Point(10000,0),
	L[tot].ag=atan2(L[tot].B.y-L[tot].A.y,L[tot].B.x-L[tot].A.x);
	L[++tot].A = Point(10000,0) , L[tot].B = Point(10000,10000),
	L[tot].ag=atan2(L[tot].B.y-L[tot].A.y,L[tot].B.x-L[tot].A.x);
	L[++tot].A = Point(10000,10000) , L[tot].B = Point(0,10000),
	L[tot].ag=atan2(L[tot].B.y-L[tot].A.y,L[tot].B.x-L[tot].A.x);
	L[++tot].A = Point(0,10000) , L[tot].B = Point(0,0),
	L[tot].ag=atan2(L[tot].B.y-L[tot].A.y,L[tot].B.x-L[tot].A.x);
	sort(L+1,L+1+tot);
	q[r++] = 1;
	for(int i=2;i<=tot;i++)
		if(dcmp(L[i].ag-L[i-1].ag))
		{
			for(;r-l>=2 && dcmp((qp[r-2]-L[i].A)*(L[i].B-L[i].A))>=0;r--);
			for(;r-l>=2 && dcmp((qp[l]-L[i].A)*(L[i].B-L[i].A))>=0;l++);
			q[r++] = i;
			qp[r-2] = Ipt(L[i].A,L[i].B,L[q[r-2]].A,L[q[r-2]].B);
		}
	for(;r-l>=2 && dcmp((qp[r-2]-L[q[l]].A)*(L[q[l]].B-L[q[l]].A))>=0;r--);
	for(;r-l>=2 && dcmp((qp[l]-L[q[r-1]].A)*(L[q[r-1]].B-L[q[r-1]].A))>=0;l++);
	if(r-l<=2) puts("0.0");
	else
	{
		qp[r-1] = Ipt(L[q[l]].A,L[q[l]].B,L[q[r-1]].A,L[q[r-1]].B);
		qp[r] = qp[l];
		double ans = 0;
		for(int i=l;i<r;i++)
		{
			ans += (qp[i]-qp[l]) * (qp[i+1]-qp[l]);
		}
		printf("%.1f\n",fabs(ans*0.5));
	}
}

2020/1/5UPD:

调整了代码长度。

20200604UPD:

$n \log n$ 半平面交在一次半平面的极角变化>=pi的时候是会有问题的，况且这个问题是不能简单修好的。