poj1279 不规则多边形的核，用半平面切割多边形 O（n^2）

最新推荐文章于 2021-05-28 18:12:57 发布

二分抄代码

最新推荐文章于 2021-05-28 18:12:57 发布

阅读量218

点赞数

分类专栏：二维点线面计算几何

本文链接：https://blog.csdn.net/liufengwei1/article/details/87372601

版权

计算几何同时被 2 个专栏收录

48 篇文章 0 订阅

订阅专栏

二维点线面

37 篇文章 0 订阅

订阅专栏

//const double inf=1e9;
//求半平面切割多边形的时候，要先放入4个极限点当整个平面，再求面积，但是由于double只能存1e300+的数，而如果inf=1e15-18就很容易在中间有乘法的计算过程爆炸损失精度，所以inf要根据题意来设置

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#define maxl 1010
#define eps 1e-8

using namespace std;

inline int sgn(double x)
{
	if(x>-eps && x<eps) return 0;
	if(x>0) return 1;
	else	return -1;
}

struct point
{
	double x,y;
	point(double a=0,double b=0)
	{
		x=a,y=b;
	}
	point operator - (const point &b)const
	{
		return point(x-b.x,y-b.y);
	}
	inline double norm()
	{
		return sqrt(x*x+y*y);
	}
};
inline double dot(const point &a,const point &b)
{
	return a.x*b.x+a.y*b.y;
}
inline double det(const point &a,const point &b)
{
	return a.x*b.y-a.y*b.x;
}
struct polygon_convex
{
	vector<point> P;
	polygon_convex(int Size=0)
	{
		P.resize(Size);
	}
};
struct halfplane
{
	double a,b,c;//ax+by+c<=0
	halfplane(point p=point(),point q=point())
	{
		a=p.y-q.y;
		b=q.x-p.x;
		c=det(p,q);
	}
	halfplane(double aa=0,double bb=0,double cc=0)
	{
		a=aa;b=bb;c=cc;
	}
};
inline double calc(halfplane &L,point &a)
{
	return a.x*L.a+a.y*L.b+L.c;
} 
point Intersect(point &a,point &b,halfplane &L)
{//use similar triangles ,line ans halfplaneline
	point res;
	double t1=calc(L,a),t2=calc(L,b);
	res.x=(t2*a.x-t1*b.x)/(t2-t1);
	res.y=(t2*a.y-t1*b.y)/(t2-t1);
	return res;
} 
polygon_convex cut(polygon_convex &a,halfplane &L)
{
	int n=a.P.size();
	polygon_convex res;
	for(int i=0;i<n;i++)
	{
		if(sgn(calc(L,a.P[i]))<0) 
			res.P.push_back(a.P[i]);
		else
		{
			int j=i-1;
			if(j<0) j=n-1;
			if(sgn(calc(L,a.P[j]))<0)
				res.P.push_back(Intersect(a.P[j],a.P[i],L));
			j=i+1;
			if(j==n) j=0;
			if(sgn(calc(L,a.P[j]))<0)
				res.P.push_back(Intersect(a.P[i],a.P[j],L));
		}	
	}
	return res;
}

const double inf=1e9;

polygon_convex core(vector<point> &a)
{
	polygon_convex res;
	res.P.push_back(point(-inf,-inf));
	res.P.push_back(point(inf,-inf));
	res.P.push_back(point(inf,inf));
	res.P.push_back(point(-inf,inf));
	int n=a.size();
	for(int i=0;i<n;i++)
	{
		halfplane L(a[i],a[(i+1)%n]);
		res=cut(res,L);
	}
	return res;
}

int n;
vector <point> a;
double ans;

inline void prework()
{
	scanf("%d",&n);
	a.resize(n);
	for(int i=0;i<n;i++)
		scanf("%lf%lf",&a[i].x,&a[i].y);
}

inline double area(vector<point> &a)
{
	double sum=0;int n=a.size();
	for(int i=0;i<n;i++)
		sum+=det(a[(i+1)%n],a[i]);
	return fabs(sum/2);
}

inline void mainwork()
{
	polygon_convex cr=core(a);
	ans=area(cr.P);
}

inline void print()
{
	printf("%.2f\n",ans);
}

int main()
{
	int t;
	scanf("%d",&t);
	for(int i=1;i<=t;i++)
	{
		prework();
		mainwork();
		print();
	}
	return 0;
}

2019.2.16 O(nlogn)更新

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#define maxl 1510
#define eps 1e-8
const double inf=1e9;
using namespace std;

inline int sgn(double x)
{
	if(x>-eps && x<eps) return 0;
	if(x>0) return 1;
	else	return -1;
}

struct point
{
	double x,y;
	point(double a=0,double b=0)
	{
		x=a,y=b;
	}
	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);
	}
	friend point operator * (const double &t,const point &a)
	{
		return point(t*a.x,t*a.y);
	}
	inline double norm()
	{
		return sqrt(x*x+y*y);
	}
};
inline double dot(const point &a,const point &b)
{
	return a.x*b.x+a.y*b.y;
}
inline double det(const point &a,const point &b)
{
	return a.x*b.y-a.y*b.x;
}
struct polygon_convex
{
	vector<point> P;
	polygon_convex(int Size=0)
	{
		P.resize(Size);
	}
};
struct halfplane
{// s->e on the left
	point s,e;
	double k;
	halfplane(point a=point(),point b=point())
	{
		s=a;e=b;
		k=atan2(e.y-s.y,e.x-s.x);
	}
	point operator &(const halfplane &b)const
    {
        double t=det(b.s-s,b.e-b.s);
		t=t/det(e-s,b.e-b.s);
		return s+t*(e-s);
    }
};
inline bool satisfy(point a,const halfplane L)
{//不允许在线上 
	return sgn(det(a-L.s,L.e-L.s))<=0;
}
inline bool HPIcmp(const halfplane &a,const halfplane &b)
{//如果平行且同向，内侧的在前 
	int res=sgn(a.k-b.k);
	return res==0 ? det(a.s-b.s,b.e-b.s)<0 : a.k<b.k;
}
halfplane Q[maxl];
void HPI(halfplane line[],int n,point res[],int &resn)
{
	int tot=n;
	sort(line,line+n,HPIcmp);
	tot=1;
	for(int i=1;i<n;i++)
	if(sgn(line[i].k-line[i-1].k)!=0)
		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(sgn(det(Q[tail].e-Q[tail].s,Q[tail-1].e-Q[tail-1].s))==0||
		   sgn(det(Q[head].e-Q[head].s,Q[head+1].e-Q[head+1].s))==0)
		   	return; 
		while(head<tail && !satisfy(Q[tail]&Q[tail-1],line[i]))
			tail--;
		while(head<tail && !satisfy(Q[head]&Q[head+1],line[i]))
			head++;
		Q[++tail]=line[i];
	}
	while(head<tail && !satisfy(Q[tail]&Q[tail-1],Q[head]))
		tail--;
	while(head<tail && !satisfy(Q[head]&Q[head+1],Q[tail]))
		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];
}


int n,resn;
point p[maxl];
halfplane a[maxl];
point res[maxl];
double ans;

inline void prework()
{
	scanf("%d",&n);
	for(int i=0;i<n;i++)
	{
		scanf("%lf%lf",&p[i].x,&p[i].y);
		if(i)
			a[i]=halfplane(p[i],p[i-1]);
	}
	a[0]=halfplane(p[0],p[n-1]);
	a[n++]=(halfplane(point(-inf,-inf),point(inf,-inf)));
	a[n++]=(halfplane(point(inf,-inf),point(inf,inf)));
	a[n++]=(halfplane(point(inf,inf),point(-inf,inf)));
	a[n++]=(halfplane(point(-inf,inf),point(-inf,-inf)));
}

inline double area(point a[],int n)
{
	double sum=0;
	if(n<3)
		return 0;
	for(int i=0;i<n;i++)
		sum+=det(a[(i+1)%n],a[i]);
	return fabs(sum/2);
}

inline void mainwork()
{
	HPI(a,n,res,resn);
	ans=area(res,resn);
}

inline void print()
{
	printf("%.2f\n",ans);
}

int main()
{
	int t;
	scanf("%d",&t);
	for(int i=1;i<=t;i++)
	{
		prework();
		mainwork();
		print();
	}
	return 0;
}