Not Only Algorithm,不仅仅是算法,关注数学、算法、数据结构、程序员笔试面试以及一切涉及计算机编程之美的内容 。。
你的位置:NoAlGo博客 » 算法 » 

最远点对问题

二维平面上有n个点,怎么寻找距离最远的两个?
前面讨论过最近点对的问题,这里研究其相对的问题,虽然题目类似,但思想方法截然不同,本文将使用graham凸包扫描法和旋转卡壳算法进行求解。

算法思路

最直观的思路是暴力枚举,枚举所有的两个点,计算距离并找出最大的。这种算法的时间复杂度为O(n^2),不够优化。

凸包(Convex Hull)是一个计算几何的概念,点集Q的凸包是指一个最小凸多边形,满足Q中的点或者在多边形边上或者在其内。直观上来讲,想象n个顶点上各有一枚钉子,使用一根橡皮圈把这所有钉子包围起来得到的图像就是这n个点的凸包。

不难想象,距离最远的点肯定在凸包的顶点上,于是可以先对所有的点求一个凸包,然后再对凸包上的点求最短距离。一般来说,凸包上的点的数量会比原始点的数量少很多,这时我们可以对凸包上的所有进行暴力枚举,所花费的时间不会太大。但对于特殊构造的数据,可能所有的n个点就构成了一个凸包,这时在凸包上枚举跟直接在原始点集上枚举是一样的,复杂度仍然为O(n^2),这时我们可以选择使用旋转卡壳方法求取最远点对。

Graham扫描法

对于凸包的形成过程,我们可以想象拿一根线,从最左下角的顶点开始,沿着逆时针方向一直往前包,最后回到最左下角这个点,线在旋转过程中碰到的顶点即组成了整个凸包。这就是Graham扫描法的整体思想。

首先需要找到最左下角的顶点,将其固定到数组的首位,然后对剩余的点进行极角排序。所谓的极角排序跟平时常见的根据横纵坐标的排序不一样,它使用每个点对于参考点的旋转角度(倾斜角)为标准进行排序,如果倾斜角一样,则距离参考点较近的点排在前面。
具体实现的时候并不适用具体的倾斜角,而是使用叉乘完成。如下图所示,A为参考点,比较B和C的顺序。只需计算向量AB和向量AC的叉乘,根据右手定则,把右手的四指从B握向C,此时大拇指朝上,表示叉乘结果为正,此时B应该排在C的前面。

    C
   /
  /
 /
A ---- B

具体的扫描过程如下图所示,从A开始,沿着逆时针进行扫描。为了在扫描过程中得到一个凸的形状,我们每次只能够向左转,如果出现向右转的时候,需要把前面的顶点删除,判断出现向左转为止。

F---------E
|        /
|       D
|        \ 
|         C
|        /
A ---- B

具体扫描过程如下:

  1. 一开始把A和B入栈。
  2. 下一点C,是AB方向的左边,C入栈。
  3. 下一点D,是BC方向的左边,D入栈。
  4. 下一点E,是CD方向的右边,不符要求,把D出栈。现在E是BC方向的左边,E入栈。
  5. 下一点F,是DE方向的左边,F入栈。
  6. 最终栈中的顶点为ABCEF,即为凸包的顶点。

具体实现过程的向左转向右转仍然是通过向量叉乘实现的,比如上图的AB和AC的叉乘为正,说明是向左转(逆时针转)。

另外还有一点是当多个点共线时,这里的做法是忽略中间的点,直接选取首尾的两个端点。具体参考代码的相关部分。

旋转卡壳算法

旋转卡壳可以用于求凸包的直径(最远点对)、宽度,两个不相交凸包间的最大距离和最小距离等。旋转卡壳是指两条相互平行的直线,如下图的上下两条线所示,沿着凸包的表面同步旋转。

----B------
   /  \
  A    D
   \  /
----C-----

考虑使用旋转卡壳求凸包的直径,即最远点对的问题,具体的步骤如下:

  1. 计算凸包Y顶点Y坐标的最大值ymax和最小值ymin。
  2. 通过ymin和ymax构造两条水平线,初始化最大距离为此时的距离。
  3. 沿着凸包旋转这两条线,直到其中一条与凸包的某条边缘重合,计算测试的距离,更新最大距离。
  4. 重复步骤3,直到再次回到ymin和ymax,此时得到的最大距离即为直径。

但直接使用上述思想进行实现的话代码非常复杂,我们使用另一种思路。对于直径上的两点A和B,做出一对平行线,然后沿着凸包表面旋转,知道其中一条和凸包的某条边重合,此时另一个顶点是距离这条直线最远的点。于是,我们可以枚举凸包上的所有边,依次求得离这条边最远的顶点,然后计算这个顶点到边两个端点的距离并取较大的一个,枚举完成时得到的最大距离即为直径长度。
直接实现时需要枚举所有边,每次又要枚举所有顶点,时间复杂度为O(n^2)。但是可以注意到,当枚举的边逆时针旋转时,最远点也是跟着逆时针变化,这样我们可以不用每次枚举所有的顶点,直接从上次的最远点开始继续计算即可,此时复杂度为O(n)。关于这个性质的严谨证明这里不做深入,具体的使用方法参考代码实现。

注意,在代码中,寻找离当前边最远的顶点时并没有直接去计算点到直线的距离,而是根据三角形面积进行判断,距离直线越远的点构成的三角形的面积越大。而三角形的面积可以通过叉乘很方便地计算出来。如下图所示,三角形ABC的面积即为AB和AC的叉乘的一半(比较时忽略的1/2的常系数)。即
S = AB x AC = (xb-xa, yb-ya) x (xc-xa, yc-ya) = (xb-xa)*(yc-ya) – (xc-xa)*(yb-ya)
其中(xa, ya)等表示点坐标,x表示叉乘,*表示数值乘法。

    C
   / \
  /   \
 /     \
A ----- B

代码实现

下面是具体的代码实现过程。

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const int mxsz = 100000 + 10; //最大的点数量
struct Point
{
	double x, y;
	Point (double a = 0, double b = 0): x(a), y(b) {}
} v[mxsz]; //顶点集
int stk[mxsz], top; //graham扫描得到的凸包顶点下标

//计算欧几里得距离
inline double calDis(const Point &p1, const Point &p2)
{
	return sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y));
}

//计算向量ab和向量cd的叉乘
double cross(const Point &a, const Point &b, const Point &c, const Point &d)
{
	double x1 = b.x - a.x, y1 = b.y - a.y, x2 = d.x - c.x, y2 = d.y - c.y;
	return x1 * y2 - x2 * y1;
}

bool cmp(const Point &a, const Point &b)
{
	double t = cross(v[0], a, v[0], b);
	if (fabs(t) < 1e-6) return calDis(v[0], a) < calDis(v[0], b);
	else return t > 0; //右手法则,叉乘为正,则由a到b为逆时针
}

//graham扫描法求凸包
void graham (int n)
{
	//最左下角的点放到首位
	int id = 0; 
	for (int i = 1; i < n; i++)
		if (v[i].y < v[id].y || (v[i].y == v[id].y && v[i].x < v[id].x))
			id = i;
	swap(v[0], v[id]);

	//极角排序
	sort(v+1, v+n, cmp);

	//前两个点直接入栈
	for (int i = 0; i <= 1; i++) stk[i] = i;
	top = 1;

	for (int i = 2; i < n; i++)
	{
		while (top > 0 && cross(v[stk[top-1]], v[stk[top]], v[stk[top-1]], v[i]) <= 0)
			top--; //共线的点只保留两端,如果n个点均共线,则最终只得到首尾两个端点
		stk[++top] = i;
	}
}

//旋转卡壳求凸包的直径
double rorating_caliper()
{
	double ans = 0;
	stk[++top] = 0; //最后一个点即为最开始的点
	for (int i = 0, j = 2; i < top; i++)
	{
		while (cross(v[stk[i]], v[stk[i+1]], v[stk[i]], v[stk[j]]) 
			 < cross(v[stk[i]], v[stk[i+1]], v[stk[i]], v[stk[j+1]]))
				j = (j + 1) % top;
		ans = max(ans, max(calDis(v[stk[i]], v[stk[j]]), calDis(v[stk[i+1]], v[stk[j]])));
	}
	return ans;
}

int main()
{
	//输入n个点
	int n; scanf("%d", &n);
	for (int i = 0; i < n; i++) scanf("%lf%lf", &v[i].x, &v[i].y);

	//graham扫描法求凸包
	graham(n);

	//旋转卡壳求最远距离
	double len = rorating_caliper();
	printf("Longest distance is: %lf\n", len);
}
上一篇: 下一篇:

我的博客

NoAlGo头像编程这件小事牵扯到太多的知识,很容易知其然而不知其所以然,但真正了不起的程序员对自己程序的每一个字节都了如指掌,要立足基础理论,努力提升自我的专业修养。

站内搜索

最新评论