对于几何相关的计算,中学的数学已经涉及了一部分。对于两个有向线段a和b,如何判断线段a在线段b的顺时针方向还是逆时针方向,两个线段是否相交。我们可以通过计算两个线段的斜率来判断它们之间的相对方向,使用方程组来求解两个线段是否有交点。这些方法计算机也可以使用,但是计算机计算很大的特点就是精度有限,计算斜率和解方程组都会有除法运算,这样的方法是不够好的。

除了中学接触的几何问题,还有线段集合的相交判断,以及寻找凸包等问题。完整代码以及测试代码位于Gist

数据结构

// point data structure
struct Point
{
	int x, y;
	Point(): Point(0, 0) {}
	Point(int x, int y): x(x), y(y) {}
	Point operator+(const Point &p2) const { return Point(x+p2.x, y+p2.y); }
	Point operator-(const Point &p2) const { return Point(x-p2.x, y-p2.y); }
	bool operator==(const Point &p2) const { return x==p2.x && y==p2.y; }
	bool operator!=(const Point &p2) const { return !((*this)==p2); }
};

#define INNER_PRODUCT(p1,p2) ((p1).x*(p2).y-(p1).y*(p2).x)

一个点也就是一个二维向量,支持加法和减法。向量的乘法有两种:一种是点乘x1x2+y1y2,另一种是叉乘x1y2-x2y1,叉乘(也就是内积)在几何计算当中非常重要。

// segment data structure
struct Segment
{
	Point start, end;
	Segment() = default;
	Segment(const Point &start, const Point &end): start(start), end(end) {}
	bool operator==(const Segment &seg2) const { return start == seg2.start && end == seg2.end; }
	bool operator!=(const Segment &seg2) const { return !((*this) == seg2); }
};

一个有向线段包含了开始顶点和结束顶点。

线段方向

内积就是两个二维向量构成的平行四边形的面积,通过内积是可以确定两条有向线段的相对方向的。根据内积的计算公式:

p1×p2=p1.xp2.y-p1.yp2.x=-(p2.xp1.y-p2.yp1.x)=-p2×p1

如果内积为正数,说明向量p1位于向量p2的顺时针方向;如果内积为负数,那么说明向量p1位于向量p2的逆时针方向;如果内积为零,那么p1和p2相向或者反向。

线段的起点通常不是原点,所以还需要提供一个共同的起点(起点不同的线段之间的比较通常没有太大意义)。

// return positive number if p1->p3 is on the clockwise direction of p1->p2, negative number if anticlockwise, zero elsewise
#define DIRECTION(p1,p2,p3)	(INNER_PRODUCT((p3)-(p1),(p2)-(p1)))

线段相交

两条线段相交

判断两个线段是否相交基于一下三种情况:

  • 如果对于两条线段中的任意一条线段,另外一条线段的起点和终点分别位于线段的两侧,可以直接确定相交。

  • 如果四个点中的某个点位于另外一条线段所在的直线上,那么通过判断这个点是否位于这条线段上来确定是否相交。

  • 不符合上述情况,不相交。

#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))
// when p and s are on the same line, return true if p on s
#define ON_SEGMENT(p,s) (MIN((s).start.x, (s).end.x) <= (p).x && (p).x <= MAX((s).start.x, (s).end.x)\
		&& MIN((s).start.y, (s).end.y) <= (p).y && (p).y <= MAX((s).start.y, (s).end.y))

bool intersect(const Segment &seg1, const Segment &seg2)
{
	int d1 = DIRECTION(seg1.start, seg1.end, seg2.start);
	int d2 = DIRECTION(seg1.start, seg1.end, seg2.end);
	int d3 = DIRECTION(seg2.start, seg2.end, seg1.start);
	int d4 = DIRECTION(seg2.start, seg2.end, seg1.end);
	if (((d1 < 0 && d2 > 0) || (d1 > 0 && d2 < 0))
		&& ((d3 < 0 && d4 > 0) || (d3 > 0 && d4 < 0)))
		return true;
	if (d1 == 0 && ON_SEGMENT(seg2.start, seg1)) return true;
	if (d2 == 0 && ON_SEGMENT(seg2.end, seg1)) return true;
	if (d3 == 0 && ON_SEGMENT(seg1.start, seg2)) return true;
	if (d4 == 0 && ON_SEGMENT(seg1.end, seg2)) return true;
	return false;
}

一组线段相交

对于一组线段,怎么判断是否存在两个以上的线段相交呢?我们需要两个数据结构:

  • 事件列表:也就是每一个线段的端点,检查相交就是按照事件进行;

  • 线段集合:需要保存当前在检查范围内的线段,按照与当前事件点所在的竖线交点的y坐标大小进行排序(也就是说,排序函数需要访问当前的事件点所在的x坐标值)。

  1. 首先将全部线段的端点加入到一个事件列表中,这个列表按照点的x坐标值排序,如果坐标值相同,左端点优先,如果都是左端点,那么按照y坐标值排序。

  2. 然后按照顺序处理事件列表中的“事件点”。

  • 如果事件点是左端点,那么将这个端点所在的线段加入到有序集合中,并且检查这个线段与前驱和后继线段是否相交;

  • 如果事件点是右端点,那么检查这个线段的前驱和后继线段是否相交,然后将这个线段从集合中删除。

接下来考虑,为什么这样是对的?考虑两条线段相交,那么可以分为两种情况:

  • 在交点的左侧,不存在其他线段处于两条线段之间,那么在第二条线段加入在集合中后,第一条线段也就是它的前驱或者后继,可以马上检测到相交。

  • 在交点的左侧,存在其他线段处于两条线段之间,那么在最后一个处于两条线段之间的线段被移除是,两条线段相交分别是被删除线段的前驱或者后继,于是就检测到了相交。

bool intersect(const std::vector<Segment> &segs)
{
	// intern class
	struct Event
	{
		bool end;
		Point point;
		Segment segment;
		Event() = default;
		Event(bool end, const Point &point, const Segment &segment): end(end), point(point), segment(segment) {}
	};
	// construct event list
	std::vector<Event> events(segs.size()*2);
	for (int i = 0; i < segs.size(); i++) {
		Point left = segs[i].start;
		Point right = segs[i].end;
		if (left.x > right.x)
			std::swap(left, right);
		events[2*i] = Event(false, left, segs[i]);
		events[2*i+1] = Event(true, right, segs[i]);
	}
	std::sort(events.begin(), events.end(), [](const Event &lhs, const Event &rhs)->bool{
		if (lhs.point.x != rhs.point.x)
			return lhs.point.x < rhs.point.x;
		if (lhs.end != lhs.end)
			return lhs.end < rhs.end;
		return lhs.point.y < rhs.point.y;
	});
	// check intersect
	int x;
	auto compare = [&x](const Segment &lhs, const Segment &rhs)->bool{
		int dx1 = lhs.end.x - lhs.start.x;
		int dx2 = rhs.end.x - rhs.start.x;
		int dy1 = lhs.end.y - lhs.start.y;
		int dy2 = rhs.end.y - rhs.start.y;
		return (dy1*(x-lhs.start.x)+dx1*lhs.start.y)*dx2 
			 < (dy2*(x-rhs.start.x)+dx2*rhs.start.y)*dx1;
	};
	std::set<Segment, decltype(compare)> segSet(compare);
	for (const Event &event : events) {
		x = event.point.x;
		if (event.end == false) {
			auto ret = segSet.insert(segSet.begin(), event.segment);
			auto precursor = std::prev(ret);
			auto successor = std::next(ret);
			if (precursor != segSet.end() && precursor != ret && intersect(*precursor, event.segment))
				return true;
			if (successor != segSet.end() && intersect(*successor, event.segment))
				return true;
		} else {
			auto ret = segSet.find(event.segment);
			auto precursor = std::prev(ret);
			auto successor = std::next(ret);
			if (precursor != segSet.end() && precursor != ret && 
				successor != segSet.end() && intersect(*precursor, *successor))
				return true;
			if (ret != segSet.end())
				segSet.erase(ret);
		}
	}
	return false;
}

寻找凸包

给定一个点的集合(点的数量肯定要大于2),要求找到一个子集,这个子集中的点构成的多边形可以围住尽可能多的点。从视觉上看,凸包可想象为一条刚好包着所有点的橡皮圈。

寻找凸包的办法有很多,下面是比较经典的两种。

Graham旋转扫除法

  1. 首先,选择点集中位于最下(如果最下方的点有多个,那么选择最左的点)的点作为开始的点p0,这个点显然是凸包中的点。

  2. 然后,将点集中剩余的按照相对于起始点p0的极角大小进行升序排序,如果两个点的大小一致,那么按照极半径大小升序排序。因为p0是最下方的点中最左的点,所以极角大小严格位于[0,Pi),可以直接使用内积确定相对次序。

  3. 最后,按照排好的顺序逐个处理剩余的点。在遍历过程中,需要用到三个点,一个最近确定为凸包中的点pi,一个待检验的点pj,以及仅次于pj的点dj+1。这个时候就会有两种情况:
    • 如果pi->pj->pj+1左转,说明直线pipj的顺时针方向以及不再有其他的点,pj确定为凸包中的点,将pj加入凸包,下一个待检验的点就是pj+1。
    • 如果如果pi->pj->pj+1没有左转,说明直线pipj的顺时针方向或者直线上有其他的点,pj不是凸包中的点,直接丢弃pj,将下一个待检验的点设为pj+1。
  4. 最后一个待检验的点一定是凸包中点,直接加入凸包。
// assert there is no duplicate points
std::vector<Point> graham(const std::vector<Point> &points)
{
	// require more than 2 points
	if (points.size() < 3) return std::vector<Point>();
	// find the downmost point (choose the leftmost in case of a tie)
	int startPointIndex = 0;
	for (int i = 1; i < points.size(); i++)
		if (points[i].y < points[startPointIndex].y 
			|| (points[i].y == points[startPointIndex].y && points[i].x < points[startPointIndex].x))
			startPointIndex = i;
	Point startPoint = points[startPointIndex];
	// sort rest points by polar angle (compare distance in case of a tie)
	std::vector<Point> restPoints(points.size()-1);
	for (int i = 0, offset = 0; i < points.size(); i++)
		if (i != startPointIndex)
			restPoints[offset++] = points[i];
	std::sort(restPoints.begin(), restPoints.end(), [&startPoint](const Point &lhs, const Point &rhs)->bool{
		int d = DIRECTION(startPoint, rhs, lhs);
		return d == 0 ? abs(lhs.x-startPoint.x) < abs(rhs.x-startPoint.x) : d > 0;
	});
	// construct convex hull
	std::vector<Point> hull = {startPoint};
	Point top = restPoints[0];
	for (int i = 1; i < restPoints.size(); i++) {
		if (DIRECTION(hull[hull.size()-1], top, restPoints[i]) < 0)	// turn left
			hull.push_back(top);
		top = restPoints[i];
	}
	hull.push_back(top);
	return hull;
}

Graham旋转扫除法中,每一个节点只会扫除一次,所以扫除过程的时间复杂度为O(n),而排序的时间复杂度O(nlgn)决定了算法最终的时间复杂度。

Jarvis步进法

在旋转扫除法当中,每次检查的目的就是让pi->pj的极角尽量小,如果有多个选择,那么pi->pj的极半径需要尽量大,其实也可以直接从起始点p0开始简单粗暴地不断寻找满足条件的点。

// assert there is no duplicate points
std::vector<Point> jarvis(const std::vector<Point> &points)
{
	// require more than 2 points
	if (points.size() < 3) return std::vector<Point>();
	// find the downmost point (choose the leftmost in case of a tie)
	Point startPoint = points[0];
	for (int i = 1; i < points.size(); i++)
		if (points[i].y < startPoint.y 
			|| (points[i].y == startPoint.y && points[i].x < startPoint.x))
			startPoint = points[i];
	// construct convex hull
	std::vector<Point> hull = {startPoint};
	Point currentPoint = startPoint;
	while (true) {
		// find the point with the minimum polar angle (choose the farthest point in case of a tie)
		bool first = true;
		Point nextPoint;
		for (const Point &point : points)
			if (point != currentPoint)
				if (first) {
					first = false;
					nextPoint = point;
				} else if (comparePolar(currentPoint, point, nextPoint)) {
					nextPoint = point;
				}
		if (nextPoint == startPoint) break;
		currentPoint = nextPoint;
		hull.push_back(currentPoint);
	}
	return hull;
}

由于算法中需要比较的极角大小位于[0,2Pi),单纯使用内积是无法比较两个有向线段的相对极角大小的,需要分情况处理。对于两个有向线段s1和s2,将两条线段的共同起点作为原点:

  • 如果s1和s2位于x轴的同一侧,那么使用内积就可以确定它们的相对极角大小,为了满足最优条件,当它们共线的时候,极半径大的线段为“小”。

  • 如果s1和s2位于x轴的两侧,当然是位于上方的线段极角比较小。

  • 如果s1和s2都和x轴重合,方向为x轴正方向的线段极角较小。如果两条线段方向相同,极半径大的线段为“小”。

  • 剩下的情况就是s1和s2有且仅有一条线段与x轴重合,说明s1和s2方向既不可能相等,它们之间的夹角角度在(0,Pi)内,所以可以直接使用内积进行比较。

// return true if the polar angle of p1->p3 is greater than p1->p2 (compare polar radius in case of a tie)
bool comparePolar(const Point &p1, const Point &p2, const Point &p3)
{
	int d2 = p2.y - p1.y;
	int d3 = p3.y - p1.y;
	if ((d2 > 0 && d3 > 0) || (d2 < 0 && d3 < 0)) {
		int d = DIRECTION(p1, p2, p3);
		return (d == 0) ? (abs(p2.x-p1.x) > abs(p3.x-p1.x)) : (d < 0);
	}
	if (d2 > 0 && d3 < 0)
		return false;
	if (d2 < 0 && d3 > 0)
		return true;
	if (d2 == 0 && d3 == 0)
		if (p2.x-p1.x > 0 && p3.x-p1.x < 0)
			return true;
		else if (p2.x-p1.x < 0 && p3.x-p1.x > 0)
			return false;	
		else return abs(p2.x-p1.x) > abs(p3.x-p1.x);	// in case of a tie
	return DIRECTION(p1, p2, p3) < 0;
}

对于每个节点,寻找下一个节点需要扫描全部的节点,如果存在h个顶点的话,时间复杂度为O(nh)。

最近节点对

如果找到点集中最近节点对?一种朴素的方法就是遍历每个节点对,时间复杂度为O(n^2)。比较经典的方法就是使用分治法:

  • 基本情况:如果一个点集大小为2或者3,那么只要简单比较就可以获得最近节点对。

  • 分解:找到一条竖线p,将点集中的点分为两部分。分别求解两个子集的最近节点对。

  • 合并:现在已经获得了两个子集的最近节点对,设两个最近节点对的最小距离为δ,如果想要存在跨越两个子集的最小节点对,那么它的距离必然小于δ。所以,首先将距离竖线超过δ的点排除在外。然后,将距离竖线δ内的点按照y值大小进行排序,接着从上到下对这些点进行检查。由于检查是从上到下进行,所以检查某个点的时候,只需要考虑在这个点下面的点。因为最小距离δ的约束,只需要检查于下面7个点构成的节点对。

在进行递归求解之前,对点集进行排序,然后在分解和合并过程中好好利用有序性质,可以让合并和分解的时间复杂度保持为O(n)。根据递归式O(T)=2*O(T/2)+O(n),算法的时间复杂度为O(nlogn)。

#define SQUARE(a) ((a)*(a))
#define DISTANCE(a,b) (SQUARE((a).x-(b).x)+SQUARE((a).y-(b).y))

std::pair<int, int> closestPair(const std::vector<Point> &points, std::vector<bool> &bitmap,
	const std::vector<int> &x, const std::vector<int> &y)
{
	if (x.size() == 2) {
		return std::make_pair(x[0], x[1]);
	} else if (x.size() == 3) {
		int mind = DISTANCE(points[x[0]], points[x[1]]);
		int d02 = DISTANCE(points[x[0]], points[x[2]]);
		int d12 = DISTANCE(points[x[1]], points[x[2]]);
		std::pair<int, int> minp = std::make_pair(x[0], x[1]);
		if (d02 < mind) {
			mind = d02;
			minp = std::make_pair(x[0], x[2]);
		}
		if (d12 < mind)
			minp = std::make_pair(x[1], x[2]);
		return minp;
	} else {
		// split x
		int mid = x.size() / 2;
		std::vector<int> lx(mid), rx(x.size()-mid);
		for (int i = 0; i < mid; i++) {
			lx[i] = x[i];
			bitmap[x[i]] = false;
		}
		for (int i = mid; i < x.size(); i++) {
			rx[i-mid] = x[i];
			bitmap[x[i]] = true;
		}
		// split y
		int lyoffset = 0, ryoffset = 0;
		std::vector<int> ly(mid), ry(y.size()-mid);
		for (int i = 0; i < y.size(); i++)
			if (bitmap[y[i]] == false)
				ly[lyoffset++] = y[i];
			else
				ry[ryoffset++] = y[i];
		// get subanswer
		std::pair<int, int> p1 = closestPair(points, bitmap, lx, ly);
		std::pair<int, int> p2 = closestPair(points, bitmap, rx, ry);
		int d1 = DISTANCE(points[p1.first], points[p1.second]);
		int d2 = DISTANCE(points[p2.first], points[p2.second]);
		int mind = d1;
		std::pair<int, int> minp = p1;
		if (d2 < d1) {
			mind = d2;
			minp = p2;
		}
		// find candidate pairs
		std::vector<int> cy;
		for (int i = 0; i < x.size(); i++)
			if (SQUARE(points[x[i]].x-points[x[mid]].x) <= mind)
				cy.push_back(x[i]);
		std::sort(cy.begin(), cy.end(), [&points](int lhs, int rhs)->bool{ return points[lhs].y < points[rhs].y; });
		for (int i = 0; i < cy.size(); i++)
			for (int j = i+1; j < cy.size() && j-i < 8; j++) {
				int d = DISTANCE(points[cy[i]], points[cy[j]]);
				if (d < mind) {
					mind = d;
					minp = std::make_pair(cy[i], cy[j]);
				}
			}
		return minp;
	}
}

std::pair<int, int> closestPair(const std::vector<Point> &points)
{
	// require more than one points
	if (points.size() < 2) throw std::invalid_argument("closestPair: require more than one points");
	// construct p, x, y, bitmap
	std::vector<int> x(points.size());
	std::vector<bool> bitmap(points.size());
	for (int i = 0; i < x.size(); i++)
		x[i] = i;
	std::vector<int> y(x);
	std::sort(x.begin(), x.end(), [&points](int lhs, int rhs)->bool{ return points[lhs].x < points[rhs].x; });
	std::sort(y.begin(), y.end(), [&points](int lhs, int rhs)->bool{ return points[lhs].y < points[rhs].y; });
	return closestPair(points, bitmap, x, y);
}