欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

计算几何_半平面求交

程序员文章站 2022-04-04 07:58:58
...

 

const double eps = 1e-8;
const double inf = 1e30;
int sign(double d){
	return d < -eps ? -1 : (d > eps);
}
struct point{
	double x, y;
	point(double _x=0, double _y=0) : x(_x), y(_y) {}
	void read(){
		scanf("%lf%lf", &x, &y);
	}
};
//直线p1-p2和直线a*x+b*y+c=0的交点
point intersect(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);
	return point((p1.x * v + p2.x * u) / (u + v), (p1.y * v + p2.y * u) / (u + v));
}
//求经过st和ed(保证这两个点不相等)的直线a*x+b*y+c=0的三个参数
void tran(point st, point ed, double& a, double& b, double& c){
	a = st.y - ed.y;
	b = ed.x - st.x;
	c = (ed.y - st.y) * st.x - (ed.x - st.x) * st.y;
}
//多边形类
struct poly{
	static const int N = 1005; //点数的最大值
	point ps[N+5]; //逆时针存储多边形的点,[0,pn-1]存储点
	int pn;  //点数
	poly() { pn = 0; }
	//加进一个点
	void push(point tp){
		ps[pn++] = tp;
	}
	//第k个位置
	int trim(int k){
		return (k+pn)%pn;
	}
	void clear(){ pn = 0; }
};

//求该凸多边形org(确保org的点是逆时针方向的)在半平面a*x+b*y+c < 0 的部分
poly incise(poly org, double a, double b, double c){
	poly ans;
	point* ps = org.ps;
	int i;
	for(i = 0; i < org.pn; i++){
		if(sign(a*ps[i].x+b*ps[i].y+c) < 1){
			ans.push(ps[i]);
		}else{
			int g = org.trim(i-1);
			if(sign(a*ps[g].x+b*ps[g].y+c) < 0){
				ans.push(intersect(ps[g], ps[i], a, b, c));
			}
			g = org.trim(i+1);
			if(sign(a*ps[g].x+b*ps[g].y+c) < 0){
				ans.push(intersect(ps[g], ps[i], a, b, c));
			}
		}
	}
	return ans;
}
//求该凸多边形org(确保org的点是逆时针方向的)在st-->ed的左边(即从st看向ed时的左手边)的部分
poly incise(poly org, point st, point ed){
	double a, b, c;
	point tp;
	tran(st, ed, a, b, c);
	//tp是st-->ed的左边的半平面的一个代表点
	tp.x = st.x - (ed.y - st.y);
	tp.y = st.y + (ed.x - st.x);
	if(sign(a*tp.x+b*tp.y+c) < 0){
		return incise(org, a, b, c);
	}else{
		return incise(org, -a, -b, -c);
	}
}
相关标签: 半平面求交