10.3 求解最近点对问题

二维空间中最近点对问题是，给定平面上的 $n$ 个点，找其中的一对点，使得在 $n$ 个点的所有点对中该点对的距离最小。这类问题在实际中有广泛的应用。例如，在空中交通控制问题中，若将飞机看作空间中移动的一个点来看待，则具有最大碰撞危险的两架飞机就是这个空间中最接近的一对点。本节介绍求解最近点对的两种算法。

10.3.1 用蛮力法求最近点对

蛮力法的过程是，分别计算每一对点之间的距离，然后找出距离最小的一对。对于给定的点集 $a$ ，采用蛮力法求 $a[leftindex,rightindex]$ 中的最近点对之间距离的算法如下：

double ClosestPoints(vector<Point> a, int leftindex, int rightindex) {
	double d, mindist = INF;
	for (int i = leftindex; i <= rightindex; ++i)
		for (int j = i + 1; j <= rightindex; ++j) {
			d = Distance(a[i], a[j]);
			if (d < mindist)
				mindist = d;
		}
	return mindist;
}

算法分析：上述算法中有两重 $for$ 循环，当求 $a[0...n-1]$ 中 $n$ 个点的最近点对时，算法的时间复杂度为 $O(n^2)$ 。

10.3.2 用分治法求最近点对

对于给定的点集 $a[0...n-1]$ ，采用分治法求最近点对距离的步骤如下：

1. 对 $a$ 中的所有点按 $x$ 坐标从小到大排序。将 $a$ 中点集复制到 $b$ 中，对 $b$ 中的所有点按 $y$ 坐标从小到大排序。设 $a$ 中最近点对距离为 $d$ 。
2. 如果 $a$ 中点数少于 $4$ ，则采用蛮力法直接计算各点的最近距离 $d$ 。
3. 求出 $a$ 中间位置的点 $a[midindex]$ ，以此位置画一条中轴线 $l$（对应的 $x$ 坐标为 $a[midindex].x$ ），将 $a$ 的所有点分割为点数大致相同的两个子集，左部分包含 $a[0...midindex]$ 的点，右部分包含 $a[midindex+1...n-1]$ 的点，同样将 $b$ 中的点相应分为两部分 $leftb$ 和 $rightb$ ，左部分称为 $S_1$（含 $a[0...minindex]$ 和 $leftb$ ，右部分称为 $S_2$（含 $a[midindex+1...n-1]$ 和 $rightb$ ），如图10.19所示。递归调用求出 $S_1$ 中点集的最近点对的距离为 $d_1$ ，递归调用求出 $S_2$ 中点集的最近点对的距离为 $d_2$ ，并求出当前最近点对的距离为 $d=\min (d_1,d_2)$ 。
4. 显然，$S_1$ 、$S_2$ 各自（独立）的任意点对之间的距离大于或等于 $d$ ，但 $S_1,S_2$ 交界的垂直带形区（由所有与中轴线的 $x$ 坐标值相差不超过 $d$ 的点构成）中的点对之间的距离可能小于 $d$ 。将 $b$ 中所有落在垂直带形区的点复制到 $b_1$ 中，对于 $b_1$ 中的任一点 $p$ ，仅需要考虑紧随 $p$ 后的 $7$ 个点，计算出从 $p$ 到这 $7$ 个点的距离，并和 $d$ 进行比较，将最小的距离存放在 $d$ 中，最后求得的 $d$ 即为 $a$ 中所有点的最近点对距离。

对于 $b_1$ 中的点 $p$ ，为什么只需要考虑紧随 $p$ 后的 $7$ 个点呢？如果 $p_L\in P_L,p_R\in P_R$ ，且 $p_L$ 和 $p_R$ 的距离小于 $d$ ，则它们必定位于以 $l$ 为中轴线的 $d\times 2d$ 的矩形区内，如图10.20所示，该区内最多有 $8$ 个点——左、右阴影正方形中最多各有 $4$ 个点，否则它们各自（内部）的最小距离就会小于 $d$ ，与 $P_L,P_R$ 各自（内部）的所有点的最小距离大于或等于 $d$ 矛盾。所以为了求得 $P_L$ 和 $P_R$ 中点之间的最小距离，只需考虑每个点 $p$ 之后的 $7$ 个点就可以了。

对于图10.13所示的点集 $a$ ，采用分治法求最近点对的过程如下：

1. 排序前的 $a[0...9]$ 为：$a_0(4,10),a_1(3,7),a_2(9,7),a_3(3,4),a_4(5,6),a_5(5,4),a_6(6,3),a_7(8,1),a_8(3,0),a_9(1,6)$ 。$a[0…9]按 $x$ 坐标排序后为：$a_9(1,6),a_1(3,7),a_3(3,4),a_8(3,0),a_0(4,10),a_4(5,6),a_5(5,4),a_6(6,3),a_7(8,1),a_2(9,7)$ 。将 $a$ 复制到 $b$ 中，$b[0...9]$ 按 $y$ 坐标排序后为：$a_8(3,0),a_7(8,1),a_6(6,3),a_3(3,4),a_5(5,4),a_9(1,6),a_4(5,6),a_1(3,7),a_2(9,7),a_0(4,10)$

2. 求出中间位置 $midindex=4$ ，对应顶点 $a_0$ ，左部分为 $(a_9,a_1,a_3,a_8,a_0)$ ，右部分为 $(a_4,a_5,a_6,a_7,a_2)$ 。中间部分为 $(a_8,a_3,a_5,a_4,a_1,a_0)$（按 $y$ 从小到大排列），其中，左部分包含中间位置的顶点。

3. 处理左部分 $(a_9,a_1,a_3,a_8,a_0)$ 。

   1. 求出中间位置 $midindex=2$ ，对应顶点 $a_3$ ，左部分为 $(a_9,a_1,a_3)$ ，右部分为 $(a_8,a_0)$ ，中间部分为 $(a_8,a_3,a_9,a_1)$ （按 $y$ 从小到大排列）。其中，左部分包含中间位置的顶点。

   2. 处理左部分 $(a_9,a_1,a_3)$ 。由于顶点个数少于 $4$ ，采用蛮力法求出最近距离 $d_{11}=2.23607$ ，对应的点对是 $a_9,a_1$ 。

   3. 处理右部分 $(a_8,a_0)$ 。由于顶点个数少于 $4$ ，采用蛮力法求出最近距离 $d_{12}=10.0499$ ，对应的点对是 $a_8,a_0$ 。
      左右部分合起来求出 $d_1=\min (d_{1}1,d_{1}2)=\min (2.23607,10.0499)=2.23607$ ，对应的点对是 $a_9,a_1$ 。

   4. 中间部分为 $(a_8,a_3,a_9,a_1)$（按 $y$ 从小到大排列）。

      • 考虑 $a_8$ ，在 $y$ 方向上后面没有小于 $d_1$ 的顶点；
      • 考虑 $a_3$ ，在 $y$ 方向上后面小于 $d_1$ 的顶点只有顶点 $a_9$ ，求出 $a_3$ 到 $a_9$ 的距离为 $2.82843$ 。
      • 考虑 $a_9$ ，在 $y$ 方向上后面小于 $d_1$ 的顶点只有顶点 $a_1$ ，求出 $a_9$ 到 $a_1$ 的距离为 $2.23607$ 。
      • 考虑 $a_1$ ，在 $y$ 方向后面没有其他顶点。

      求出中间部分的最近距离 $d_{13}=2.23607$ ，对应的点对是 $a_9,a_1$ 。

   这样合并得到左部分 $(a_9,a_1,a_3,a_8,a_0)$ 的结果 $d_1=\min (d_1,d_{13})=(2.23607,2.23607)=2.23607$ ，对应的点对是 $a_9$ 和 $a_1$ 。

4. 处理右部分 $(a_4,a_5,a_6,a_7,a_2)$ 。

   1. 求出中间位置 $midindex=7$ ，对应顶点 $a_6$ ，左部分为 $(a_4,a_5,a_6)$ ，右部分为 $(a_7,a_2)$ ，中间部分为 $(a_6, a_5, a_4)（
      按 $y$ 从小到大排列）。其中，左部分包含中间位置的顶点。

   2. 处理左部分 $(a_4,a_5,a_6)$ 。由于顶点个数少于 $4$ ，采用蛮力法求出最近距离 $d_{21}=1.41421$ ，对应的点对是 $a_5,a_6$ 。

   3. 处理右部分 $(a_7,a_2)$ 。由于顶点个数少于 $4$ ，采用蛮力法求出最近距离 $d_{22}=6.08276$ ，对应的点对是 $a_7,a_2$ 。
      左右部分求出 $d_2=\min (d_{21},d_{22})=\min (1.41421,6.08276)=1.41421$ ，对应的点对是 $a_5,a_6$ 。

   4. 中间部分为 $(a_6,a_5,a_4)$（按 $y$ 从小到大排列）。

      • 考虑 $a_6$ ，在 $y$ 方向后面小于 $d_2$ 的顶点只有顶点 $a_5$ ，求出 $a_6$ 到 $a_5$ 的距离为 $1.41421$ 。
      • 考虑 $a_5$ ，在 $y$ 方向后面没有小于 $d_2$ 的顶点。
      • 考虑 $a_4$ ，在 $y$ 方向后面没有其他顶点。

      求出中间部分的最近距离 $d_{23}=1.41421$ ，对应的点对是 $a_6,a_5$ 。

   这样合并得到右部分 $(a_4,a_5,a_6,a_7,a_2)$ 的结果 $d_2=\min (d_2,d_{23})=(1.41421,1.41421)=1.41421$ ，对应的点对是 $a_5,a_6$ 。

5. 考虑左右部分，求出 $d=\min (d_1,d_2)=1.41421$ 。中间部分点集为 $(a_8,a_3,a_5,a_4,a_1,a_0)$（按 $y$ 从小到大排列）。

   1. 考虑 $a_8$ ，在 $y$ 方向上后面没有小于 $d$ 的顶点。
   2. 考虑 $a_3$ ，在 $y$ 方向上后面小于 $d$ 的顶点只有顶点 $a_5$ ，求出 $a_3$ 到 $a_5$ 的距离为 $2$ 。
   3. 考虑 $a_5$ ，在 $y$ 方向上后面没有小于 $d$ 的顶点。
   4. 考虑 $a_4$ ，在 $y$ 方向上后面小于 $d$ 的顶点只有顶点 $a_1$ ，求出 $a_4$ 到 $a_1$ 的距离为 $2.23607$ 。
   5. 考虑 $a_1$ ，在 $y$ 方向上后面没有小于 $d$ 的顶点。
   6. 考虑 $a_0$ ，在 $y$ 方向上后面没有其他顶点。

   求出中间部分 $(a_8,a_3,a_5,a_4,a_1,a_0)$ 的最近距离 $d_3=2$ ，对应的点对是 $a_3,a_5$ 。

6. 合并最终结果：$d=\min (d,d_3)=\min (1.41421,2)=1.41421$ ，对应的点对是 $a_5,a_6$ 。
   

采用分治法求最近点对的算法如下：

bool pointxcmp(Point &p1, Point &p2) { // 用于点按x坐标递增排序
    return p1.x < p2.x;
}
bool pointycmp(Point &p1, Point &p2) { // 用于点按y坐标递增排序
    return p1.y < p2.y;
}
double ClosestPoints2(vector<Point> &a, vector<Point> b, int leftindex, int rightindex) { 
    // 递归求a[leftindex...rightindex]中最近点对的距离
    vector<Point> leftb, rightb, b1;
    int i, j, midindex;
    double d1, d2, d3 = INF, d;
    if ((rightindex - leftindex + 1) <= 3) {	// 少于4个点,直接用蛮力法求解
    	d = ClosestPoints(a, leftindex, rightindex);
		return d;
    }
    midindex = (leftindex + rightindex) / 2;	// 求中间位置
    for (i = 0; i < b.size(); ++i) {			// 将b中点集分为左右两部分
    	if (b[i].x < a[midindex].x) 
	   		leftb.push_back(b[i]);				// leftb中所有点的x<中轴x,且按y值从小到大排序
		else
	   		rightb.push_back(b[i]);				// rightb中所有点的x>=中轴x,且按y值从小到大排序
    }
    d1 = ClosestPoints2(a, leftb, leftindex, midindex);
    d2 = ClosestPoints2(a, rightb, midindex + 1, rightindex);
    d = min(d1, d2);							// 当前最小距离d=min(d1,d2)
    // 求中间部分点对的最小距离
    for (i = 0; i < b.size(); ++i) {			// 将b中间宽度为2d的带状区域内的子集,复制到b1中
    	if (fabs(b[i].x - a[midindex].x <= d))
	  		b1.push_back(b[i]);					// b1中的点按y值从小到大排序
	}
    double tmpd3;
    for (i = 0; i < b1.size(); ++i) {			// 求b1中最近点对
    	for (j = i + 1; j < b1.size(); ++j) {	// 对于点b1[i],这样的点b1[j]最多7个
		    if (b1[j].y - b1[i].y >= d) break;
		    tmpd3 = Distance(b1[i], b1[j]);
		    if (tmpd3 < d3)
		    	d3 = tmpd3;
		}
    }
    d = min(d, d3);
    return d;
}
double ClosestPoints1(vector<Point> &a) {
    // 求a中最近点对的距离
    printf("排序前:\n");
    for (Point &v : a) v.disp();
    printf("\n");

    sort(a.begin(), a.end(), pointxcmp);     	// 按x坐标从小到大排序
    printf("按 x 坐标排序后:\n");
    for (Point &v : a) v.disp();
    printf("\n");
    
    vector<Point> b(a.begin(), a.end());     	// 将a中点集复制到b中去
    sort(b.begin(), b.end(), pointycmp);     	// 按y坐标从小到大排序    
    printf("按 y 坐标排序后:\n");
    for (Point &v : a) v.disp();
    printf("\n");
    return ClosestPoints2(a, b, 0, a.size() - 1);
}

算法分析：在求 $a[0...n-1]$ 中 $n$ 个点的最近点时，设执行时间为 $T(n)$ ，求左、右部分中最近点对的时间为 $T(n/2)$ ，求中间部分的时间为 $O(n)$ ，则：

T(n) = O(1)					当n < 4时 
T(n) = 2T(n / 2) + O(n) 	其他情况

从而推出算法的时间复杂度为 $O(n{\log }_{2}n)$ 。