计算几何学 | 线段相交问题 | 平面扫描 | Segment Intersections: Manhattan Geometry | C/C++实现 | 曼哈顿几何

问题描述

现给出n条平行于x轴或y轴的线段，请输出其交点数。

输入：
第1行输入线段数n。接下来n行输入n条线段。每条线段按照下述格式给出：
$x_1$ $y_1$ $x_2$ $y_2$
上面分别为线段两端点的坐标。各输入数据均为整数。
输出：
输出交点的总数，占1行。
限制：
1 ≤ n ≤ 100000
互相平行的2条或更多线段之间不存在重叠的点或线段
交点数不超过1000000。

输入示例

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

输出示例

3

讲解

求n条线段的交点，可以用抽选配对的方法遍历所有可能的线段对，将交点一一列举（计数）。但是复杂度为$O(n^2)$，不理想。

与轴平行的线段相交问题（曼哈顿几何）可以通过平面扫描（sweep）高效地求解。平面扫描算法的思路是将一条与x轴（或y轴）平行的直线向上（向右）平行移动，在移动过程中寻找交点。这条直线称为扫描线。

扫描线并不是按照固定的间隔逐行扫描，而是在每次遇到平面上线段的端点时停止移动，然后检查该位置上的线段交点。为了进行上述处理，我们需要先将输入的线段端点按照y值排序，让扫描线向y轴正方向移动。

在扫描线移动的过程中，算法会将扫描线穿过的垂直线段（与y轴平行）临时记录下来，等到扫描线与水平线段（与x轴平行）重叠时，检查水平线段的范围内是否存在垂直线段上的点，然后将这些点作为交点输出。为提高处理效率，我们可以应用二叉搜索树来保存扫描穿过的垂直线段。线段相交问题的平面扫描算法具体如下

平面扫描：
1.将已输入线段的端点按y坐标升序排列，添加至表EP
2.将二叉搜索树T置为空
3.按顺序取出EP的端点（相当于让扫描线自下而上移动），进行以下处理：
如果取出的端点为垂直线段的上端点，则从T中删除该线段的x坐标
如果取出的端点为垂直线段的下端点，则将该线段的x坐标插入T
如果取出的端点为水平线段的左端点（扫描线与水平线段重合时），将该水平线段的两端点作为搜索范围，输出T中包含的值（即垂直线段的x坐标）。

使用平衡的二叉搜索树后，1此搜索的复杂度为$O(logn)$，由于这个值小于2n，所以二叉树带来的复杂度为$O(nlogn)$。算法整体的复杂度还与交点数k有关，因此本题中介绍的平面扫描算法的复杂度为$O(nlogn+k)$。

另外这个问题还可用区间树（segment tree）高效求解。

平面搜索：

#define BOTTOM 0
#define LEFT 1
#define RIGHT 2
#define TOP 3

class EndPoint {
	public:
		Point p;
		int seg, st;//输入线段的ID，端点的种类
		EndPoint() {}
		EndPoint(Point p, int seg, int st): p(p), seg(seg), st(st) {}
		
		bool operator < (const EndPoint &ep) const {
			//按y坐标升序排序
			if(p.y == ep.p.y) {
				return st < ep.st;//y相同时，按照下端点、左端点、右端点、上端点的顺序排列 
			} else return p.y < ep.p.y; 
		}
}; 

EndPoint EP[2 * 100000];//端点列表

//线段相交问题：曼哈顿几何
int manhattanIntersection(vector<Segment> S) {
	int n = S.size();
	
	for(int i = 0, k = 0; i < n; i++) {
		//调整端点p1、p2，保证左小右大
		if(S[i].p1.y == S[i].p2.y) {
			if(S[i].p1.x > S[i].p2.x) swap(S[i].p1, S[i].p2);
		} else if (S[i].p1.y > S[i].p2.y) swap(S[i].p1, S[i].p2);
		
		if(S[i].p1.y == S[i].p2.y) {
			EP[k++] = EndPoint(S[i].p1, i, LEFT);
			EP[k++] = EndPoint(S[i].p2, i, RIGHT);
		} else {
			EP[k++] = EndPoint(S[i].p1, i, BOTTOM);
			EP[k++] = EndPoint(S[i].p2, i, TOP);
		}
	}
	
	sort(EP, EP + (2 * n));//按端点的y坐标升序排列
	
	set<int> BT;//二叉搜索树
	BT.insert(1000000001);//设置标记
	int cnt = 0;
	
	for(int i = 0; i < 2 * n; i++) {
		if(EP[i].st == TOP) {
			BT.erase(EP[i].p.x);//删除上端点 
		} else if(EP[i].st == BOTTOM) {
			BT.insert(EP[i].p.x);//添加下端点 
		} else if(EP[i].st == LEFT) {
			set<int>::iterator b = BT.lower_bound(S[EP[i].seg].p1.x);//O(log n)
			set<int>::iterator e = BT.upper_bound(S[EP[i].seg].p2.x);//O(log n)
			cnt += distance(b, e);//加上b和e的距离（点数） O(k) 
		}
	} 
	
	return cnt;
} 

实现上述平面扫描算法时要注意各种处理的顺序，以免在一条扫描线上同时进行多个处理时遗漏交点。上述算法在一条扫描线上同时进行线段（x坐标的值）的删除、插入、搜索时，会按照下端点、左端点、右端点、上端点的顺序排列端点，从而避免这一问题。

AC代码如下

#include<iostream>
#include<cmath>
#include<vector> 
#include<algorithm>
#include<set>
using namespace std;

#define EPS (1e-10)
#define equals(a, b) (fabs((a) - (b)) < EPS)

class Point {//Point类，点 
	public:
		double x, y;
		
		Point(double x = 0, double y = 0): x(x), y(y) {}

		Point operator + (Point p) { return Point(x + p.x, y + p.y); }
		Point operator - (Point p) { return Point(x - p.x, y - p.y); }
		Point operator * (double a) { return Point(a * x, a * y); }
		Point operator / (double a) { return Point(x / a, y / a); }

		double abs() { return sqrt(norm()); }
		
		double norm() { return x * x + y * y; }
		
		bool operator < (const Point &p) const {
			return x != p.x ? x < p.x : y < p.y;
		}

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

typedef Point Vector;//Vector类，向量 

struct Segment{//Segment 线段 
	Point p1, p2;
};

double dot(Vector a, Vector b) {//内积 
	return a.x * b.x + a.y * b.y;
}

#define BOTTOM 0
#define LEFT 1
#define RIGHT 2
#define TOP 3

class EndPoint {
	public:
		Point p;
		int seg, st;//输入线段的ID，端点的种类
		EndPoint() {}
		EndPoint(Point p, int seg, int st): p(p), seg(seg), st(st) {}
		
		bool operator < (const EndPoint &ep) const {
			//按y坐标升序排序
			if(p.y == ep.p.y) {
				return st < ep.st;//y相同时，按照下端点、左端点、右端点、上端点的顺序排列 
			} else return p.y < ep.p.y; 
		}
}; 

EndPoint EP[2 * 100000];//端点列表

//线段相交问题：曼哈顿几何
int manhattanIntersection(vector<Segment> S) {
	int n = S.size();
	
	for(int i = 0, k = 0; i < n; i++) {
		//调整端点p1、p2，保证左小右大
		if(S[i].p1.y == S[i].p2.y) {
			if(S[i].p1.x > S[i].p2.x) swap(S[i].p1, S[i].p2);
		} else if (S[i].p1.y > S[i].p2.y) swap(S[i].p1, S[i].p2);
		
		if(S[i].p1.y == S[i].p2.y) {
			EP[k++] = EndPoint(S[i].p1, i, LEFT);
			EP[k++] = EndPoint(S[i].p2, i, RIGHT);
		} else {
			EP[k++] = EndPoint(S[i].p1, i, BOTTOM);
			EP[k++] = EndPoint(S[i].p2, i, TOP);
		}
	}
	
	sort(EP, EP + (2 * n));//按端点的y坐标升序排列
	
	set<int> BT;//二叉搜索树
	BT.insert(1000000001);//设置标记
	int cnt = 0;
	
	for(int i = 0; i < 2 * n; i++) {
		if(EP[i].st == TOP) {
			BT.erase(EP[i].p.x);//删除上端点 
		} else if(EP[i].st == BOTTOM) {
			BT.insert(EP[i].p.x);//添加下端点 
		} else if(EP[i].st == LEFT) {
			set<int>::iterator b = BT.lower_bound(S[EP[i].seg].p1.x);//O(log n)
			set<int>::iterator e = BT.upper_bound(S[EP[i].seg].p2.x);//O(log n)
			cnt += distance(b, e);//加上b和e的距离（点数） O(k) 
		}
	} 
	
	return cnt;
} 

int main(){
	int n;
	cin>>n;
	
	vector<Segment> S;
	Segment s;
	
	while(n--){
		cin>>s.p1.x>>s.p1.y>>s.p2.x>>s.p2.y;
		S.push_back(s);
	}
	
	cout<<manhattanIntersection(S)<<endl;
}

注：以上本文未涉及代码的详细解释参见：计算几何学