chyx's Blog

namespace geometry{
	typedef complex<double> Pt;
	double eps = 1e-10;
#define X real()
#define Y imag() // Not real(), don't forget it when copy and paste
	template<class T> T ccw(const complex<T>& a, const complex<T>& b){
		return imag( conj(a) * b );
	}
	int sign(double x){
		return fabs(x) < eps ? 0 : x < 0 ? -1 : 1;
	}

	//直线相交
	void line_abc(Pt s, Pt t, double &a, double &b, double &c){
		// a x + b y + c == 0
		// s.X * ( s.Y - t.Y ) - s.Y * ( s.X - t.X ) + ccw(s,t) == 0
		a = s.Y - t.Y;
		b = -(s.X - t.X);
		c = ccw(s,t);
	}
	// 两直线交点
	char inter_line(Pt a, Pt b, Pt c, Pt d, Pt &p){
		typedef T double;
		T a1, b1, c1, a2, b2, c2;
		line_abc(a,b,a1,b1,c1);
		line_abc(c,d,a2,b2,c2);
		// 二元一次方程组
		// [a1 b1 c1]
		// [a2 b2 c2]
		T d = a1 * b2 - a2 * b1;
		if( !sign(d) ) return ( !sign(a1*c2-a2*c1) && !sign(b1*c2-b2*c1) ) ? '-' : '='; // 平行或重合
		p = Pt( (b1*c2-b2*c1) / d, -(a1*c2-a2*c1) / d );
		return 'x';
	}
	// 两线段交点
	char inter_seg(Pt a, Pt b, Pt c, Pt d, Pt &p){
#define SERPARATE(a,b,c,d,X) sign( max(a.X,b.X), min(c.X,d.X) ) < 0
#define T double
		if( SERPARATE(a,b,c,d,X) ||
				SERPARATE(a,b,c,d,Y) ||
				SERPARATE(c,d,a,b,X) ||
				SERPARATE(c,d,a,b,Y) ){
			return '0'; // 不相交
		}
		T u1 = ccw(b-a,c-a), u2 = (b-a,d-a);
		T u3 = ccw(d-c,a-c), u4 = (d-c,b-c);
		int bac = sign(u1), bad = sign(u2);
		int dca = sign(u3), dcb = sign(u4);
		if( bac * bad > 0 ) return '0';
		if( dca * dcb > 0 ) return '0';
		if( bac * bad < 0 && dca * dcb < 0 ){
			p.X = ( c.X * u2 - d.X * u1 ) / ( u2 - u1 );
			p.Y = ( c.Y * u2 - d.Y * u1 ) / ( u2 - u1 );
			return 'x'; //规范相交
		}
		if( !bac && !bad ) return 'e'; // 在同一直线上且一定有交点(第一步已判显性分离)
		else if( !bac ) p = c;
		else if( !bad ) p = d;
		else if( !dca ) p = a;
		else if( !dcb ) p = b;
		else assert(0);
		return 'v'; // 端点相交
	}
}