泰森多边形半平面求交 - 洛谷

欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。

往期相关背景半平面求交点击前往
voronoi 图求解点击前往

题目大意

题目链接
https://www.luogu.com.cn/problem/P3297

小杨家是一个矩阵，所有亲戚都在矩阵里。小杨一开始只有离一个亲戚最近。
小杨在家里可以随意移动，每次移动都会被最近的亲戚监视到，问小杨想要离开家，最少要被几个亲戚监视。

解析

整体思路是先对所有亲戚所在位置进行三角剖分，然后求出泰森多边形。

这样可以得到每个点所管辖区域邻近区域的关系。

那么从一个区域走到另一个区域是增加1个人头。

最后利用广度优先搜索最短路即可求解。

那么核心问题就是如何求得每个点与附近点的关系。

在这里插入图片描述

最直观的做法是，先进行三角剖分，由于voronoi区域是通过所有边的垂直平分线包围形成的。所以三角形的边就代表点之间的连通关系。

但是会有一些特殊情况，

在这里插入图片描述
如上图，P1, P2, P3, P4是共圆的，蓝色线是其三角剖分。

本来红色线是一个分割线，由于红色线与P3P2的交点与4点距离一样，所以如果想直接从P3走到P2那么会被4个人同时监视到，所以该分割失效。

改进算法：
通过计算垂直平分线，得到区域边界。

利用半平面求交，依次加入边（同时建立对应点到点的关系），如果之前算出来的交点在加入的边上，需要删除之前的边，同时之前建立点的关系也要删除。

注意: 半平面求交，对于非封闭区域需要从第1条边开始往队列里加入。

特殊情况：没有人看管，少于3人看管，看管人员共线

代码


#include<stdio.h>
#include<cmath>
#include <algorithm>
#include <vector>
#include <list>
#include <cstring>
#include <utility>
#include <queue>
#include <map>


using namespace std;
const double EPS = 1e-8;

const int N = 1e6 + 10;
const int M = 1e6 + 10;


int cmp(double d) {
	if (abs(d) < EPS)return 0;
	if (d > 0)return 1;
	return -1;
}

class Point {
public:
	double x, y;
	int id;

	Point() {}
	Point(double a, double b) :x(a), y(b) {}
	Point(const Point& p) :x(p.x), y(p.y), id(p.id) {}

	void in() {
		scanf("%lf %lf", &x, &y);
	}
	void out() {
		printf("%f %f\n", x, y);
	}

	double dis() {
		return sqrt(x * x + y * y);
	}

	double dis2() {
		return x * x + y * y;
	}

	Point operator -() const {
		return Point(-x, -y);
	}

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

	Point operator +(const Point& p) const {
		return Point(x + p.x, y + p.y);
	}
	Point operator *(double d)const {
		return Point(x * d, y * d);
	}

	Point operator /(double d)const {
		return Point(x / d, y / d);
	}


	void operator -=(Point& p) {
		x -= p.x;
		y -= p.y;
	}

	void operator +=(Point& p) {
		x += p.x;
		y += p.y;
	}
	void operator *=(double d) {
		x *= d;
		y *= d;
	}

	void operator /=(double d) {
		this ->operator*= (1 / d);
	}

	bool operator<(const Point& a) const {
		return x < a.x || (abs(x - a.x) < EPS && y < a.y);
	}

	bool operator==(const Point& a) const {
		return abs(x - a.x) < EPS && abs(y - a.y) < EPS;
	}
};

// 向量操作

double cross(const Point& a, const Point& b) {
	return a.x * b.y - a.y * b.x;
}

double dot(const Point& a, const Point& b) {
	return a.x * b.x + a.y * b.y;
}


class Point3D {
public:
	double x, y, z;

	Point3D() {}
	Point3D(double a, double b, double c) :x(a), y(b), z(c) {}
	Point3D(const Point3D& p) :x(p.x), y(p.y), z(p.z) {}

	double dis() {
		return sqrt(x * x + y * y + z * z);
	}

	double dis2() {
		return x * x + y * y + z * z;
	}

	Point3D operator -(const Point3D& p) const {
		return Point3D(x - p.x, y - p.y, z - p.z);
	}

	Point3D operator +(const Point3D& p) const {
		return Point3D(x + p.x, y + p.y, z + p.z);
	}
	Point3D operator *(double d)const {
		return Point3D(x * d, y * d, z * d);
	}

	Point3D operator /(double d)const {
		return Point3D(x / d, y / d, z / d);
	}


	void operator -=(Point3D& p) {
		x -= p.x;
		y -= p.y;
		z -= p.z;
	}

	void operator +=(Point3D& p) {
		x += p.x;
		y += p.y;
		z += p.z;
	}
	void operator *=(double d) {
		x *= d;
		y *= d;
		z *= d;
	}

	void operator /=(double d) {
		this ->operator*= (1 / d);
	}
};

// 向量操作
Point3D cross(const Point3D& a, const Point3D& b) {
	return Point3D(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x,
		a.x * b.y - a.y * b.x);
}

double dot(const Point3D& a, const Point3D& b) {
	return a.x * b.x + a.y * b.y + a.z * b.z;
}


class Line {
public:
	Point front, tail;
	int ind;
	double ang;
	Line() {}
	Line(const Point& a, const Point& b) :front(a), tail(b) {
		ang = atan2(front.y - tail.y, front.x - tail.x);
	}
	Line(const Point& a, const Point& b, int i) :front(a), tail(b), ind(i) {
		ang = atan2(front.y - tail.y, front.x - tail.x);
	}

};

/*
0 不相交
1 相交
0 平行/重合
*/
int cross(const Line& a, const Line& b) {
	Point dir1 = a.front - a.tail;
	Point dir2 = b.front - b.tail;
	if (cmp(cross(dir1, dir2)) == 0) {
		return 0;
	}

	if (cmp(cross(a.front - b.tail, dir2)) * cmp(cross(a.tail - b.tail, dir2)) >= 0)return 0;
	if (cmp(cross(b.front - a.tail, dir1)) * cmp(cross(b.tail - a.tail, dir1)) >= 0)return 0;
	return 1;
}


int inCircle(Point p0, Point p1, Point p2, Point p3) {
	Point d1 = p1 - p0;
	Point d2 = p2 - p0;
	if (cross(d1, d2) < 0)return inCircle(p0, p2, p1, p3); // 保证平面法向向上

	// 构建映射点
	Point3D lift0(p0.x, p0.y, p0.dis2());
	Point3D lift1(p1.x, p1.y, p1.dis2());
	Point3D lift2(p2.x, p2.y, p2.dis2());
	Point3D lift3(p3.x, p3.y, p3.dis2());

	Point3D z1(lift1 - lift0), z2(lift2 - lift0);
	Point3D normal = cross(z1, z2); // 计算平面法向
	double project = dot(normal, lift3 - lift0); // 计算点到平面距离

	return cmp(project);
}



class EdgeDelaunay {
public:
	int id;
	std::list<EdgeDelaunay>::iterator c;
	EdgeDelaunay(int id = 0) { this->id = id; }
};

class Delaunay {
public:
	std::list<EdgeDelaunay> head[N];  // graph
	Point p[N];
	int n = 0;

	void init(int psize, Point ps[]) {
		this->n = psize;
		for (int i = 0; i < psize; ++i) head[i].clear();
		memcpy(this->p, ps, sizeof(Point) * n);
		std::sort(this->p, this->p + n);
		divide(0, n - 1);
	}

	void addEdge(int u, int v) {
		head[u].push_front(EdgeDelaunay(v));
		head[v].push_front(EdgeDelaunay(u));
		head[u].begin()->c = head[v].begin();
		head[v].begin()->c = head[u].begin();
	}

	void divide(int l, int r) {
		if (r - l <= 1) {  // #point <= 2
			for (int i = l; i <= r; i++)
				for (int j = i + 1; j <= r; j++) addEdge(i, j);
			return;
		}
		int mid = (l + r) / 2;
		divide(l, mid);
		divide(mid + 1, r);

		std::list<EdgeDelaunay>::iterator it;
		int nowl = l, nowr = r;

		for (int update = 1; update;) {
			// 查找左边最低线位置
			update = 0;
			Point ptL = p[nowl], ptR = p[nowr];
			for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
				Point t = p[it->id];
				double v = cross(ptL - ptR, t - ptR);
				if (cmp(v) > 0 || (cmp(v) == 0 && (t - ptR).dis() < (ptL - ptR).dis())) {
					nowl = it->id, update = 1;
					break;
				}
			}
			if (update) continue;
			// 查找右边最低线位置
			for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
				Point t = p[it->id];
				double v = cross(ptR - ptL, t - ptL);
				if (cmp(v) < 0 || (cmp(v) == 0 && (t - ptL).dis() < (ptL - ptR).dis())) {
					nowr = it->id, update = 1;
					break;
				}
			}
		}

		addEdge(nowl, nowr);  // 添加基线

		for (; true;) {
			Point ptL = p[nowl], ptR = p[nowr];
			int ch = -1, side = 0;
			for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
				if (cmp(cross(ptR - ptL, p[it->id] - ptL)) <= 0)continue; // 判断夹角是否小于180
				if (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0) {
					ch = it->id, side = -1;
				}
			}
			for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
				if (cmp(cross(p[it->id] - ptR, ptL - ptR)) <= 0) continue;// 判断夹角是否小于180
				if (ch == -1 || inCircle(ptL, ptR, p[ch], p[it->id]) < 0) {
					ch = it->id, side = 1;
				}
			}
			if (ch == -1) break;  // 所有线已经加完
			if (side == -1) {
				for (it = head[nowl].begin(); it != head[nowl].end();) {
					// 判断是否相交，边缘不算相交
					if (cross(Line(ptL, p[it->id]), Line(ptR, p[ch]))) {
						head[it->id].erase(it->c);
						head[nowl].erase(it++);
					}
					else {
						it++;
					}
				}
				nowl = ch;
				addEdge(nowl, nowr);
			}
			else {
				for (it = head[nowr].begin(); it != head[nowr].end();) {
					// 判断是否相交，边缘不算相交
					if (cross(Line(ptR, p[it->id]), Line(ptL, p[ch]))) {
						head[it->id].erase(it->c);
						head[nowr].erase(it++);
					}
					else {
						it++;
					}
				}
				nowr = ch;
				addEdge(nowl, nowr);
			}
		}
	}

	std::vector<std::pair<int, int> > getEdge() {
		std::vector<std::pair<int, int> > ret;
		ret.reserve(n);
		std::list<EdgeDelaunay>::iterator it;
		for (int i = 0; i < n; i++) {
			for (it = head[i].begin(); it != head[i].end(); it++) {
				ret.push_back(std::make_pair(p[i].id, p[it->id].id));
			}
		}
		return ret;
	}
};


int cmp(const Line& a, const Line& b) {
	//auto ta = atan2(a.front.y - a.tail.y, a.front.x - a.tail.x);
	//auto tb = atan2(b.front.y - b.tail.y, b.front.x - b.tail.x);

	return cmp(a.ang - b.ang);
}


// 点在直线哪一边>0 左边，<0边
int SideJudge(const Line& a, const Point& b) {
	return cmp(cross(a.front - a.tail, b - a.tail));
}


int LineSort(const Line& a, const Line& b) {
	int c = cmp(a, b);
	if (c)return c < 0;
	return SideJudge(b, a.front) > 0;
}

/*
点p 到 p+r 表示线段1
点q 到 q+s 表示线段2
线段1 上1点用 p' = p+t*r (0<=t<=1)
线段2 上1点用 q' = q+u*s (0<=u<=1)
让两式相等求交点 p+t*r = q+u*s
两边都叉乘s
(p+t*r)Xs = (q+u*s)Xs
pXs + t*rXs = qXs
t = (q-p)Xs/(rXs)
同理，
u = (p-q)Xr/(sXr) -> u = (q-p)Xr/(rXs)

以下分4种情况：
1. 共线，sXr==0 && (q-p)Xr==0, 计算 (q-p)在r上的投影在r长度上的占比t0，
计算(q+s-p)在r上的投影在r长度上的占比t1，查看[t0, t1]是否与范围[0，1]有交集。
如果t0>t1, 则比较[t1, t0]是否与范围[0，1]有交集。
t0 = (q-p)*r/(r*r)
t1 = (q+s-p)*r/(r*r) = t0 + s · r / (r · r)
2. 平行sXr==0 && (q-p)Xr!=0
3. 0<=u<=1 && 0<=t<=1 有交点
4. 其他u, t不在0到范围内，没有交点。
*/
pair<double, double> intersection(const Point& q, const Point& s, const Point& p, const Point& r) {
	// 计算 (q-p)Xr
	auto qpr = cross(q - p, r);
	auto qps = cross(q - p, s);

	auto rXs = cross(r, s);
	if (cmp(rXs) == 0)return { -1, -1 }; // 平行或共线
	// 求解t, u
	// t = (q-p)Xs/(rXs)
	auto t = qps / rXs;

	// u = (q-p)Xr/(rXs)
	auto u = qpr / rXs;

	return { u, t };
}

Point LineCross(const Line& a, const Line& b) {
	Point dira = a.front - a.tail;
	Point dirb = b.front - b.tail;

	auto p = intersection(a.tail, dira, b.tail, dirb);
	return a.tail + dira * p.first;
}


class HalfPlane {
public:
	vector<Line> lines;

	void addLine(const Line& a) {
		lines.push_back(a);
	}

	vector<int> run(vector<Point>& ps) {
		sort(lines.begin(), lines.end(), LineSort);
		vector<int> q(lines.size() + 10);
		vector<Point> t(lines.size() + 10);

		// 查找到第1条线
		int k = lines.size() - 1;
		bool isBoundary = false;
		for (; k >= 0; --k) {
			int nk = (k + 1)%lines.size();
			if (cmp(cross(lines[k].front - lines[k].tail, lines[nk].front - lines[nk].tail)) <= 0) {
				isBoundary = true;
				k = nk;
				break;
			}
		}

		if (!isBoundary)k = 0;

		int l = -1, r = 0;
		q[0] = k;
		for (int i = (k+1)%lines.size(); i !=k; i=(i+1)%lines.size()) {
			while (r - l > 1 && SideJudge(lines[i], t[r]) <= 0)r--;
			while (r - l > 1 && SideJudge(lines[i], t[l + 2]) <= 0)l++;
			q[++r] = i;
			t[r] = LineCross(lines[q[r]], lines[q[r - 1]]);
		}

		while (r - l > 2 && SideJudge(lines[q[l + 1]], t[r]) <= 0)r--;
		vector<int> ans;
		for (int i = l + 1; i <= r; ++i) {
			ans.push_back(lines[q[i]].ind);
			if(i>l+1)ps.push_back(t[i]);
		}
		
		if (!isBoundary) {
			t[r + 1] = LineCross(lines[q[l + 1]], lines[q[r]]);
			ps.push_back(t[r + 1]);
		}

		return ans;
	}
};


class Edge {
public:
	Edge() {}
	Edge(int t, int n, double w) :to(t), next(n), weight(w) {}
	int index;
	int from;
	int to;
	int next;
	double weight;
	bool isValid() const { return to >= 0; }
};

class Graph {
public:
	int size;
	Graph() {}
	void init(int n) {
		size = n;
		head.assign(n, -1);
		edge.resize(M * 2);
		len = 0;
		emptyEdge = Edge(-1, -1, 0);
	}

	Edge headEdge(int a) {
		if (head[a] < 0)return emptyEdge;
		return edge[head[a]];
	}

	Edge nextEdge(const Edge& ed) {
		if (ed.next < 0)return emptyEdge;
		return edge[ed.next];
	}

	void add(int a, int b, double w) {
		//printf("add : %d %d %.6lf\n", a, b, w);
		edge[len] = Edge(b, head[a], w);
		edge[len].index = len;
		head[a] = len;
		len++;
	}

	void del(int ind) {
		edge[ind].to = -1;
	}

private:
	vector<int> head;
	int len;
	Edge emptyEdge;
	vector<Edge> edge;
};

Graph g;

Point oiPs[N];
Delaunay de;
bool isBoundary[N];

int bfs(int n, int start) {
	map<int, bool> vis;
	int step = 1;

	queue<int> qu;
	qu.push(start);
	vis[start] = true;

	while (!qu.empty()) {
		int len = qu.size();
		while (len--) {
			auto front = qu.front();
			qu.pop();
			if (isBoundary[front]) return step;

			for (Edge ed = g.headEdge(front); ed.isValid(); ed = g.nextEdge(ed)) {
				if (vis[ed.to])continue;
				vis[ed.to] = true;
				qu.push(ed.to);
			}
		}

		step++;
	}

	return step;
}

void  solve() {
	int n, t;
	scanf("%d", &t);
	int x1, y1, x0, y0;
	while (t--) {
		scanf("%d", &n);
		scanf("%d%d%d%d", &x1, &y1, &x0, &y0);
		// 输入亲戚坐标
		for (int i = 0; i < n; ++i) {
			int a, b;
			scanf("%d%d", &a, &b);
			oiPs[i].x = a;
			oiPs[i].y = b;
			oiPs[i].id = i;
		}
		if (n == 0) {
			puts("0");
			continue;
		}
		if (n < 4) {
			puts("1");
			continue;
		}

		// 所有点共线
		bool oneLine = true;
		for (int i = 2; i < n; ++i) {
			if (cmp(cross(oiPs[0] - oiPs[1], oiPs[0] - oiPs[i])))oneLine = false;
		}
		if (oneLine) {
			puts("1");
			continue;
		}

		// if (n > 600) exit(1);
		de.init(n, oiPs);

		auto oiedges = de.getEdge();
		vector<vector<int>> link(n, vector<int>());
		for (auto oie : oiedges) {
			link[oie.first].push_back(oie.second);
		}

		g.init(n);

		int start = 0;
		for (int i = 0; i < n; ++i) {
			// 遍历每个边界点，收集邻边，并按照逆时针排序。
			int len = link[i].size();
			auto &ind = link[i];
			HalfPlane hp;
			// 求voronoi 边界之间交点
			for (int j = 0; j < len; ++j) {
				Point mid = (oiPs[i] + oiPs[ind[j]]) / 2;
				Point dir = oiPs[ind[j]] - oiPs[i];
				dir = { -dir.y, dir.x }; // 旋转90度

				hp.addLine(Line(mid+dir, mid, j));
			}
			vector<Point> ps;
			auto ids = hp.run(ps);
			sort(ids.begin(), ids.end());
			int k = 0;
			for (int j = 0; j < len && k<ids.size(); ++j) {
				if (ids[k] == j) {
					g.add(i, ind[j],1);

					k++;
				}
			}

			isBoundary[i] = ps.size()!=ids.size();
			for (auto& p : ps) {
				if (isBoundary[i])break;
				//if (p.x > 0 && p.x < x1 && p.y>0 && p.y < y1)continue;
				if(p.x<0 || p.x>x1 || p.y<0 ||p.y>y1) isBoundary[i] = true;
			}
			if ((oiPs[i] - Point(x0, y0)).dis() < (oiPs[start] - Point(x0, y0)).dis())start = i;
		}
		
		printf("%d\n", bfs(n, start));

	}
}



int main() {
	solve();
	return 0;

}


/*
2
4
10 10 5 5
5 6
3 5
7 5
5 3
17
14 12 7 6
7 11
6 9
7 7
1 10
2 20
1 6
2 6
1 1
2 2
5 1
5 2
13 1
12 2
12 7
13 7
12 11
13 11

5
10 10 5 5

3 3
5 5
3 7
7 3
7 7


5
10 10 5 5

0 0
5 5
0 10
10 0
10 10

3 
10 10 5 5

0 0
5 5
0 10

4
10 10 5 5

0 0
5 5
10 10
6 6



10
10 10 5 5

5 5
4 4
6 4
4 6
6 6
7 7
7 3
3 7
7 3
1 1


*/