计算几何_半平面求交

 

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);
	}
}