计算几何模板(部分)

点的定义

struct Point {
	double x, y;
	Point(double x = 0, double y = 0) : x(x), y(y) {}
	bool operator < (const Point& rhs) const {
		return x < rhs.x;
	}
};

Point operator + (const Point& A, const Point& B) {
	return Point(A.x + B.x, A.y + B.y);
}

Point operator - (const Point& A, const Point& B) {
	return Point(A.x - B.x, A.y - B.y);
}

Point operator * (double k, const Point& A) {
	return Point(A.x * k, A.y * k);
}

ostream& operator << (ostream& out, const Point& p) {
	out << "(" << p.x << "," << p.y << ")";
	return out;
}

//点积
double Dot(const Point& p1, const Point& p2) {
	return p1.x * p2.x + p1.y * p2.y;
}
//向量的模长(长度)
double Length(Point p) {
	return sqrt(Dot(p, p));
}
//求两条线段p0p1到p0p2的夹角
int Angle(const Point& p0, const Point& p1, const Point& p2) {
	double d = Dot((p1 - p0), (p2 - p0));
	if (d == 0) return 0; //直角
	else if (d > 0) return 1; //锐角
	else return -1; //钝角
}
//矢量叉积
double Det(const Point& p1, const Point& p2) {
	return p1.x * p2.y - p2.x * p1.y;
}
//通过叉积判断两直线之间的相对方向 
int Direction(const Point& p0, const Point& p1, const Point& p2) {
	double d = Det((p1 - p0), (p2 - p0));
	if (d == 0) return 0; // 三点共线
	else if (d > 0) return 1; //p1p0在p2p0的顺时针方向
	else return -1; //p1p0在p2p0的逆时针方向
}

//两点之间的距离
double Distance(const Point& A, const Point& B) {
	return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
//点到线段的距离(不是直线) p0到线段p1p2的距离
double DistPtoSegment(Point p0, Point p1, Point p2) {
	Point v1 = p2 - p1, v2 = p1 - p2, v3 = p0 - p1, v4 = p0 - p2;
	if (p1.x == p2.x && p1.y == p2.y) { //两点重合变成算两点之间的距离
		return Length(p0 - p1);
	} 
	if (Dot(v1, v3) < 0) {
		return Length(v3);
	}
 	else if (Dot(v2, v4) < 0) {
		return Length(v4);
	}else {
		return fabs(Det(v1, v3)) / Length(v1);
	}
}
//判断点p0是否在p1和p2表示的矩形内(p1p2为矩形对角线)
bool InRectangle(const Point& p0, const Point& p1, const Point& p2) {
	return Dot((p1 - p0), (p2 - p0)) <= 0;
}
//判断点p0是否在线段p1p2上 p0在直线P1P2上且在以p1p2为对角线的矩形内
bool OnSegment(Point p0, Point p1, Point p2) {
	return Dot(p1 - p0, p2 - p0) && Dot((p1 - p0), (p2 - p0)) <= 0;
}
//判断两条线段是否平行, 叉积为0即平行(共线)
bool Parallel(Point p1, Point p2, Point p3, Point p4) {
	return Det((p2 - p1), (p4 - p3)) == 0;
}

//判断线段p1p2和p3p4是否相交,只要点p1和p2两个点在线段p3p4的两边且点p3和p4在线段p1p2的两边，那么两条线段必然相交
bool SegIntersect(Point p1, Point p2, Point p3, Point p4) {
	int d1 = Direction(p3, p1, p4);
	int d2 = Direction(p3, p2, p4);
	int d3 = Direction(p1, p3, p2);
	int d4 = Direction(p1, p4, p2);
	if (d1 * d2 < 0 && d3 * d4 < 0) return true;
	if (d1 == 0 && OnSegment(p1, p3, p4)) return true;
	else if (d2 == 0 && OnSegment(p2, p3, p4)) return true;
	else if (d3 == 0 && OnSegment(p3, p1, p2)) return true;
	else if (d4 == 0 && OnSegment(p4, p1, p2)) return true;
	else return false;
}

struct Line {
	Point p1, p2;
	Line(Point p1, Point p2) : p1(p1), p2(p2) {}
	Line() {};
};
//线段是否相交
bool SegIntersect(const Line& l1, const Line& l2) {
	return SegIntersect(l1.p1, l1.p2, l2.p1, l2.p2);
}

//判断点是否在多边形内
//从目标点出发引一条射线，看这条射线和多边形所有边的交点数目。如果有奇数个交点，则说明在内部，如果有偶数个交点，则说明在外部。
bool PointInPolygon(Point p0, Point* a, int n) {
//判断点p0是否在点集a所形成的多边形内
	int cnt = 0;
	double x;
	for (unsigned i = 0; i < n; ++i) {
		Point p1 = a[i % n];
		Point p2 = a[(i + 1) % n]; //取多边形的一条边p1p2
		if (OnSegment(p0, p1, p2)) return true;
		if (p1.y == p2.y) continue; //水平线跳过
		if (p0.y < p1.y && p0.y < p2.y) continue; //p0在p1p2下方,跳过
		if (p0.y >= p1.y && p0.y >= p2.y) continue; //在上方,跳过
		//求交点坐标x的值
		x = (p0.y - p1.y) * (p2.x - p1.x) / (p2.y - p1.y) + p1.x; 
		if (x > p0.x) cnt++;
	}
	return (cnt % 2 == 1);
}
//求三个点构成的三角形的面积
//S(p1, p2, p3) = 以p1 - p2和p3 - p2为向量构成的平行四边形面积的一半
//=(p1 - p2) x (p3 - p2) / 2 记得取绝对值
double TriangleArea(Point p1, Point p2, Point p3) {
	return fabs(Det(p1 - p2, p3 - p2)) / 2;
}

//直线ax + by + c = 0
struct Line {
	double a, b, c;
};

//线段类型
struct Segment {
	Point s, e;
};
//求过p1, p2的直线方程
// a * p1.x + b * p1.y + c = 0
//过p1(x1, y1), p2(x2, y2)的直线方程为
//(y2 - y1)x + (x1 - x2)y + (x2y1 - x1y2) = 0
Line LineForm(const Point& p1, const Point& p2) {
	Line l;
	l.a = p2.y - p1.y;
	l.b = p1.x - p2.x;
	l.c = p2.x * p1.y - p1.x * p2.y;
	return l;
}

//计算直线方程l1和l2相交的交点
Point LineIntersect(const Line& l1, const Line& l2) {
	Point p;
	double d = l1.a * l2.b - l2.a*l1.b;
	p.x = (l2.c * l1.b - l1.c * l2.b) / d;
	p.y = (l2.a * l1.c - l1.a * l2.c) / d;
	return p;
}

//凸包Graham算法
Point p0; // 起点,最下且偏左的点
bool pointcmp(const Point& a, const Point& b) {
	int k = Direction(p0, a, b);
	if (k > 0) {
		return true;
	} else if (k == 0 && Distance(p0, a) < Distance(p0, b)) {
		return true;
	}
	return false;
}
//葛立恒算法求凸包O(nlogn)
int Graham(Point* a, int n, Point ch[]) {
	int top = -1, k = 0; //栈顶指针-1表示空栈,k表示最下且偏左的点的下标
	for (int i = 1; i < n; ++i) { //找最下且偏左的点a[k]
		if (a[i].y < a[k].y || (a[i].y == a[k].y && a[i].x < a[k].x)) {
			k = i;
		}
	}
	swap(a[0], a[k]); p0 = a[0]; //将起点a[0]放入p0中
	sort(a, a + n, pointcmp); //按极角从小到大排序,极角相同则离起点近的在前面
	ch[++top] = a[0];
	ch[++top] = a[1];
	ch[++top] = a[2]; //前三个点先进栈,因为凸包最少三个点
	for (int i = 3; i < n; ++i) { //判断其余所有点
		while(top >= 0 && Direction(ch[top - 1], a[i], ch[top]) > 0 || 
		Direction(ch[top - 1], a[i], ch[top]) == 0 && Distance(ch[top - 1], a[i]) > Distance(ch[top - 1], ch[top])) {
			top--; //存在右拐关系,栈顶元素出栈
		}
		ch[++top] = a[i]; //不存在右拐关系则进栈
	}
	return top + 1;
}
//求解平面最远点对问题或最近点对问题
//1.暴力法O(n^2)
double Mostdistp(vector<Point> a, int &maxindex1, int &maxindex2) {
	double d, maxdist = 0.0;
	for (int i = 0; i < a.size(); ++i) {
		for (int j = i + 1; j < a.size(); ++j) {
			d = Distance(a[i], a[j]);
			if (d > maxdist) {
				maxdist = d;
				maxindex1 = i;
				maxindex2 = j;
			}
		}
	}
	return maxdist;
}
//旋转卡壳法求最远点对
//旋转卡壳法的基本思想是,对于给定的点集,先采用Graham扫描法求出来一个凸包a
//然后根据凸包上的每一条边找到离他最远的一个点
//亦即卡着外壳旋转一圈,这便是旋转卡壳法的由来
//最远点对肯定是位于凸包上的某两个点
double RotatingCalipersl(Point ch[], int m, int &maxindex1, int &maxindex2) {
	double maxdist = 0.0;
	int j = 1; //距离某个处理边最远的点在ch数组中的下标
	for (int i = 0; i < m; ++i) {
		while(fabs(Det(ch[i] - ch[(i + 1) % m], ch[(j + 1) % m] - ch[(i + 1) % m])) > 
			  fabs(Det(ch[i] - ch[(i + 1) % m], ch[j] - ch[(i + 1) % m]))) {
			j = (j + 1) % m; //以面积来判断,面积大说明要远一些
		}
			double d1 = Distance(ch[i], ch[j]);
			if (d1 >maxdist) {
				maxdist = d1;
				maxindex1 = i;
				maxindex2 = j;
			}
			double d2 = Distance(ch[(i + 1) % m], ch[j]);
			if (d2 > maxdist) {
				maxdist = d2;
				maxindex1 = i + 1;
				maxindex2 = j;
			}
	}
	return maxdist;
}
//旋转卡壳算法O(nlogn)
double RotatingCalipers(Point* a, int n, Point ch[], Point& p1, Point& p2) {
	int index1 = 0, index2 = 0;
	int m = Graham(a, n, ch);
	double maxdist = RotatingCalipersl(ch, m, index1, index2);
	p1 = ch[index1];
	p2 = ch[index2];
	return maxdist;
}
//n个点找三个点构成一个三角形,三角形最大面积, n个点的凸包中点的个数
double RotatingCalipers(Point ch[], int m, int &index1, int &index2, int &index3) {
	double maxArea = 0.0;
	int i = 0;
	int j = i + 1; 
	int k = j + 1;
	do{ //遍历ch数组,先固定一个点ch[i]
		int p = i, q = j, r = k;
		while (fabs(Det(ch[i] - ch[j], ch[i] - ch[k])) < fabs(Det(ch[i] - ch[j], ch[i] - ch[(k + 1) % m]))) {
			k = (k + 1) % m;
		}
		if (maxArea < fabs(Det(ch[i] - ch[j], ch[i] - ch[k]))) {
			maxArea = fabs(Det(ch[i] - ch[j], ch[i] - ch[k]));
			index1 = i;
			index2 = j;
			index3 = k;
		}
		while (fabs(Det(ch[i] - ch[j], ch[i] - ch[k])) < fabs(Det(ch[i] - ch[(j + 1) % m], ch[i] - ch[k]))) {
			j = (j + 1) % m;
		}
		if (maxArea < fabs(Det(ch[i] - ch[j], ch[i] - ch[k]))) {
			maxArea = fabs(Det(ch[i] - ch[j], ch[i] - ch[k]));
			index1 = i;
			index2 = j;
			index3 = k;
		}
		while (fabs(Det(ch[i] - ch[j], ch[i] - ch[k])) < fabs(Det(ch[(i + 1) % m] - ch[j], ch[(i + 1) % m] - ch[k]))) {
			i = (i + 1) % m;
		}
		if (maxArea < fabs(Det(ch[i] - ch[j], ch[i] - ch[k]))) {
			maxArea = fabs(Det(ch[i] - ch[j], ch[i] - ch[k]));
			index1 = i;
			index2 = j;
			index3 = k;
		}
		if (i == p && j == q && k == r) {
			k = (k + 1) % m;
		}
	}while(k > 0);
	return maxArea * 0.5;
}