最邻近点对问题(Closest-Pair Problem)：二维的分治解法详解

这是该系列的第二篇。一维、三维的最邻近点对问题（Closest-Pair Problem）：戳这里
PS：建议先快速浏览一维问题的分治解法。因为各维度的解决思路具有高度的关联性。

---

1 问题描述

现在二维下，点不再是线性分布了。对于两个二维点$p_i=(x_i,y_i)$和$p_j=(x_j,y_j)$，它们之间的欧几里得距离为$\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}$。

暴力算法只能枚举每个不同的点对。它的时间复杂度是$O(n^2)$。

但是我们可以沿用一维的思维模式来解决这个问题。并且这时，由于排序也需要$O(nlogn)$，按照分治思想可以把暴力法的复杂度降到$O(nlogn)$。

2 算法描述

2.1 开始前

首先对于一个乱序的点集$P$，我们需要分别对$x$和$y$进行排序，得到$P_x$和$P_y$（注意是分别，而不是先后）。前者用于divide，后者用于merge。

2.2 Divide步骤

▲ 第一步：拆分

对于任意的点集$P$，其中点的数量$|P|=n$，类似于一维情况，可以$P_x$中第$\lfloor n/2\rfloor$个点的坐标（记为$x^{\ast }$）处，将其分为$Q$区与$R$区。

准确来说，分区的标准是点在$P_x$中的位置，小于等于$\lfloor n/2\rfloor$划入$Q$区，反之$R$区。不过按$x^{\ast }$来理解会比较方便。

▲ 第二步: 维护有序数组

对于$Q$区和$R$区，我们想要维护4个有序数组：按$x$或$y$排好序的$Q_x$， $Q_y$， $R_x$，$R_y$。其中，按$x$排序的数组，$Q_x$和$R_x$，按索引就可以分开了。

而$Q_y$和$R_y$，需要新建数组，顺次遍历$P_y$并判断从属的区域来构建，花费$O(n)$时间，且它们自然就是有序的。注意，如果你想通过$Q_x$和$R_x$重排序得到，会影响整个算法的复杂度。

现在原问题就变成了两个完全相同的子问题。最终只有小于等于3个点时，可以直接比较每个点对得出答案。

2.3 Merge步骤

假设$Q$区找到了局部最邻近点 { q 1 , q 2 } \{q_1, q_2\} {q1​,q2​}，距离为${\delta }_1$，同理$R$区找到了 { r 1 , r 2 } \{r_1, r_2\} {r1​,r2​}，距离为${\delta }_2$。并且在所有被分割线割离的点对中，即所有点对 { p i , p j } \{p_{i}, p_{j}\} {pi​,pj​} 且 $p_i\in Q\wedge p_j\in R$，有一最小距离${\delta }_3$。现在对于原问题的解，有完全相同于一维的3种情况：

• 最小距离的点对在被分割线割离的点对集合中，即${\delta }_3$是最小，应返回对应被割离的点对。
• $Q$区中找到的点对距离最小，应返回 { q 1 , q 2 } \{q_1, q_2\} {q1​,q2​}。
• $R$区中找到的点对距离最小，应返回 { r 1 , r 2 } \{r_1, r_2\} {r1​,r2​}。

现在算法的关键点就在于找出${\delta }_3$。一维情况被割离的点对中，需要比较的只有一个，而二维情况则大不一样，需要枚举两个区域的点的所有组合来比较，每次Merge时都将消耗$O(n^2)$的时间。

▲ 能否利用${\delta }_1$与${\delta }_2$，减少寻找${\delta }_3$的开销呢？

毕竟，我们只需要找可能比$Q$和$R$中局部最近更近的点对。如果令$\delta =min({\delta }_1,{\delta }_2)$，把$x^{\ast }$左右宽为$\delta$的区域记为$S$区，那么我们只需要检查这个区域内的点即可。

遗憾的是，这样做并没有减少算法的开销上限。 如果所有点都集中在这个区域里，那么仍然不可避免检查所有点对。这时，就要用到我们一直维护的$P_y$了。

假设我们要检查点$s$，它属于$Q$区。以它为起点，取出$S$区域中高为$\delta$的一块区域。然后如上图所示，我们可以把这块$\delta \times 2\delta$的区域划分为8个box。一个box内最远距离是对角线距离$\delta /\sqrt{2}<\delta$，这意味着每个box中至多存在一点。那么根据鸽巢原理（Pigeonhole Principle），对于点$s$，仅需要检查它与比它稍大的7个点之间的距离即可。超出7个点后，距离一定会超过这块区域的高度$\delta$，因而可以不用考虑。（为什么仅检查稍大的部分？原因是遍历时，稍小的部分已经检查过了）

▲ 具体做法：

从小到大遍历有序的$P_y$，如果点在距$x^{\ast }$小于$\delta$的区域内，则将它加入$S$，这样$S$也是有序的。遍历$S$，每次只比较当前点和在其之上的7个点之间的距离。这个过程需要$7n=O(n)$的时间。

PS：如果分别构造了$S_1$和$S_2$（三维中必须要这么做），那么还需要额外记录每一个点在另一半区域中的参考位置。这样，就可以仅比较对应的半个区域内的4个点来寻找${\delta }_3$。

3 算法分析

Divide步骤的拆分和维护是$O(n)$的，Merge步骤也是$O(n)$的，易知消耗的递推公式为：

$$f_n=\{\begin{array}{rlr}& 2f_{n/2}+O(n)& n>3\\ & O(1)& n\le 3\end{array}$$

由主定理可知总消耗为$O(nlogn)$。排序消耗$O(nlogn)$，则整个算法也是$O(nlogn)$。

4 伪代码

这里贴上伪码，重在理解。源码请移步 文末链接 下一节。

5 Java实现

（算了，我发现手机上我自己都不会去点连接看源码，还是放这里了好了，长就长点）
_(:з」∠)_

Point2.java

package util;

import java.util.Objects;

public class Point2 {
	public int idx;
	public long x, y;

	public Point2(long x, long y) {
		this.x = x;
		this.y = y;
	}

	/**
	 * Check whether the first point is smaller in lexicographical order
	 */
	public static boolean smaller(Point2 p1, Point2 p2) {
		if (p1.x == p2.x) {
			return p1.y < p2.y;
		}
		return p1.x < p2.x;
	}

	/**
	 * Get distance of two points
	 */
	public static double getDis(Point2 p1, Point2 p2) {
		long tmp1 = p1.x - p2.x;
		long tmp2 = p1.y - p2.y;
		return Math.sqrt(tmp1 * tmp1 + tmp2 * tmp2);
	}

	@Override
	public String toString() {
		return "Point{" +
				"idx=" + idx +
				", x=" + x +
				", y=" + y +
				'}';
	}

	@Override
	public boolean equals(Object o) {
		if (this == o) return true;
		if (o == null || getClass() != o.getClass()) return false;
		Point2 point2 = (Point2) o;
		return x == point2.x &&
				y == point2.y;
	}

	@Override
	public int hashCode() {
		return Objects.hash(x, y);
	}
}

Dim2.java

package solve;

import util.Point2;

import java.util.Arrays;
import java.util.Comparator;

import static util.Point2.getDis;

public class Dim2 {
	private Point2[] px, py;
	private Point2[] sy;    // save sy to reuse it in all sub procedure

	/**
	 * Solve the problem
	 * @param points the point set
	 * @return the point pair with smallest distance
	 */
	public Point2[] solve(Point2[] points) {
		int n = points.length;

		/* before start searching */

		px = new Point2[n];
		py = new Point2[n];
		System.arraycopy(points, 0, px, 0, n);
		System.arraycopy(points, 0, py, 0, n);
		Arrays.sort(px, Comparator.comparingLong(o -> o.x));
		Arrays.sort(py, Comparator.comparingLong(o -> o.y));

		sy = new Point2[n]; // the point in S, ranked by y

		// record the point index in px
		for (int i = 0; i < n; i++) {
			px[i].idx = i;
		}

		/* Search the ans */

		Point2[] res = find(0, n - 1, py);

		/* output the result by lexicographical order */

		if (Point2.smaller(res[0], res[1])) {
			return new Point2[]{res[0], res[1]};
		} else {
			return new Point2[]{res[1], res[0]};
		}
	}

	/**
	 * Find the pair with smallest distance
	 * @param x1 the left bound of px to find
	 * @param x2 the right bound of px to find (include)
	 * @param py the arr px[x1:x2] ranked by y
	 * @return the pair expressed in point index
	 */
	private Point2[] find(int x1, int x2, Point2[] py){
		switch (x2 - x1 + 1) {
			case 2:
				return new Point2[]{px[x1], px[x2]};
			case 3:
				double dis12 = getDis(px[x1], px[x1 + 1]);
				double dis23 = getDis(px[x1 + 1], px[x2]);
				double dis13 = getDis(px[x1], px[x2]);
				if (dis12 < dis23) {
					if (dis12 < dis13) {
						return new Point2[]{px[x1], px[x1 + 1]};
					} else {
						return new Point2[]{px[x1], px[x2]};
					}
				} else {
					if (dis23 < dis13) {
						return new Point2[]{px[x1 + 1], px[x2]};
					} else {
						return new Point2[]{px[x1], px[x2]};
					}
				}
		}

		/* Generate Qx, Rx, Qy, Ry */

		int mi = (x1 + x2) / 2;
		int idx1 = 0;
		int idx2 = 0;
		Point2[] qy = new Point2[mi - x1 + 1];
		Point2[] ry = new Point2[x2 - mi];

		for (Point2 p : py) {
			if (p.idx <= mi) {
				qy[idx1++] = p;
			} else {
				ry[idx2++] = p;
			}
		}

		/* Search recursively */

		Point2[] left = find(x1, mi, qy);
		Point2[] right = find(mi + 1, x2, ry);

		double dis1 = getDis(left[0], left[1]);
		double dis3 = getDis(right[0], right[1]);
		double delta = Math.min(dis1, dis3);

		/* Find minimum distance in crossing-area pair */

		// Generate S
		long x = px[mi].x;
		int cnt = 0;
		for (Point2 p : py) {
			if (x - delta <= p.x && p.x <= x + delta) {
				sy[cnt++] = p;
			}
		}

		Point2[] pairMin = new Point2[2];
		double disMin = delta;

		for (int i = 0; i < cnt; i++) {
			Point2 syi = sy[i];
			Point2 syj;

			// check one side
			for (int j = i + 1; j < cnt && j <= i + 7; j++) {
				syj = sy[j];
				double tmp = getDis(syi, syj);
				if (disMin > tmp) {
					disMin = tmp;
					pairMin[0] = syi;
					pairMin[1] = syj;
				}
			}
		}

		/* Compare and return */

		if (pairMin[0] != null) {
			return pairMin;
		} else if (dis1 < dis3) {
			return left;
		} else {
			return right;
		}
	}
}

下一篇：最邻近点对问题(Closest-Pair Problem)：三维的分治解法详解

---

写在最后：

1. 该问题系列的所有实现代码、暴力代码以及对拍/测试代码(Java)：戳这里
2. 本人初来乍到，这是第一次在CSDN分享博文，如有谬误欢迎指出。