平面点对问题（分治法的经典问题）

1. 问题介绍

平面上有n个点points[n]（points[i] = (x_i, y_i)），求它们之间的最近距离。

2. 常规解法

首先想到利用握手法则，两两求点之间的距离，然后进行比较可以找到最近距离，时间复杂度为O(n²)。
代码如下：

//        使用暴力法验证一下
        long t2 = System.currentTimeMillis();
        dis = Point.distance(p[0],p[1]);
        for (int i = 0; i < p.length - 1; i++) {
            for (int j = i+1; j < p.length; j++) {
                dis = Math.min(dis, Point.distance(p[j], p[i]));
            }
        }
        diff = System.currentTimeMillis() - t2;
        System.out.println(dis + "用时：" + diff + " ms");

在150000个点时，用时为43473ms，40多秒。

3. 分治法

分治法的做法是：

1. 先将点对按照x值排序；

2. 然后按x值从中间将n个点一分为二，分别递归计算左右两边的最近距离，然后经过比较求出最小值delta；

3. delta就是我们要求的值吗，不一定，因为最近的距离可能是一个点在左半区域P，另一个点在右半区域Q，所以我们还要再求这一种情况下的解。如果是按照对左边的每一个点，都将其与右边的所有点计算距离，那么复杂度公式为W(n) = 2W(n/2) + f(n)，其中f(n)表示对分别处在左右两个半区间里的点的处理，这里左边点为n/2,右边也是n/2，因此f(n) = (n/2)² = n²/4，即时间复杂度递推方程为W(n) = 2W(n/2) + n²/4，当n = 3时递归可以直接终止W(3) = 2,W(2) = 1.最终我们根据递推树得到时间复杂度仍然为O(n²)和暴力法一样。从图中可以看到可能的改进思路在于减小f(n)的规模。
   
   因此我们要想办法降低左右两边点的求解次数。注意到，因为我们已经求出了左右两个区域内部点的最小距离delta，所以对于中间线两边delta范围之外的点就不用再考虑了，因为delta范围外的点在计算跨区间距离的时候，肯定是大于delta的。
   那么对于区域内的点，以左半区域P为例，对于P内任一点p(x_p,y_p)，在右半区域Q中，我们只需要计算y坐标值在[y_p-delta, y_p+delta]之内的点即可，原因也是一样的，如果超过了这个范围，那么再与p点求的距离肯定大于delta。并且可以知道这个范围内最多只能存在6个点。见下图：
   
   我们只需要计算右边蓝色虚线与中间线围起来的区域内的点即可。其实如果是纯粹的数学计算，我们应该是以p为圆心，δ为半径画圆，只需要计算在圆内且处于右半边的点与p的距离即可。但是这是计算机，为了实现起来简单一些，就使用的矩形作为划分。下面来证明这个宽为δ，高为2δ的矩形区域内最多只能存在6个点：

   1. 如图所示大矩形被划分为了6个小矩形，小矩形的对角线很容易得出长度为：5δ/6，这个值小于δ；
   2. 假如这个大矩形内超过6个点，则必有至少一个小矩形内有2个或更多点，然而由1可知，一个矩形内的所有点间的最大距离都是小于δ的，而又因为δ是这个区域内的最小值，不可能有比δ更近的点对存在，因此出现了矛盾。所以最多只存在6个点。

   这样的话，对于左边界区的每一个点p，只需要最多在右边界区找6个点，计算距离。
   怎么找这6个点呢，最坏的情况是右边的点没有顺序，则需要线性扫描一遍，而在最坏的情况下，所有的点都在边界区内，那么仍然需要对于每一个点都计算一遍，和刚才的O(n²)是一样的。
   如果是将右边的点按y坐标值进行排序的话，排序需要O（$$\frac{n}{2}log\frac{n}{2}$$），找这6个点就快多了，仅需log(n/2)（二分查找）。这个时候的时间复杂度为$$\frac{n}{4}log\frac{n}{4}logn$$,也即O($$nlo{g}^{2}n$$)：见下图
   

4. 算法性能分析

看起来并没有减少多少复杂度？其实不是的，当n很大时，比如n为2的31次方，大概是21亿(2E9)，暴力法的数量级为:$$10^{18}$$
分治法：$$10^9\times 31\times 31$$
很显然差了7个数量级，很吓人了。所以说logn真是个好算法。

5. 进一步改进

可以看到，制约分治法性能的主要原因在于每次递归处理都需要对点按照y进行排序，如果能够提前排好序，再进行递归，那么时间复杂度可以进一步降低到$$nlogn$$。
尽管快排已经是最好的排序方法了，但是它和分治法是不兼容的，还有一种也是nlogn级别的排序方法，也是利用分治法实现的，那就是归并排序。那么我们可以将归并排序引入到这里来：

1. 递归计算左边的最小距离；
2. 递归计算右边的最小距离；
3. 左右两边合并处理，同时对左右两边按照y值进行归并排序，这一步需要的时间为O(n)的量级。但是需要额外的辅助空间，相当于用空间换取时间。

或者还有一种方法是根据屈婉玲教授的视频课上讲的那样，预处理降低f(n)的复杂度。具体做法就是：

1. 使用两个数组X，Y，X按x值排序，Y按照y值排序；
2. 在划分时，X就直接从中间划开；
3. 对于Y，则需要遍历一遍，对于其中的每个点，如果按照x值划分到左边，则就划到左边，否则就划到右边，即Y数组的划分也是按照x值来的。这样做，划分后，保证了两边按y值仍然是有序的。排序只需要在最开始进行一次，递归过程中不再需要，在递归过程中的每次划分，其时间消耗为O(n)。

这两种改进的方法都可以做到在进行合并处理的时候将时间复杂度降到O(n)级别，最终的时间复杂度为：

6. 算法实现

代码有点多，但是里面有一些是我为了练习java，写着玩的。java刚学了不到一个月，嘻嘻嘻。

import java.util.*;

public class FindNearstPointPair {
    public static void main(String[] args) {
        Point p[] = new Point[1500000];
        Random ran = new Random(1);
        for (int i = 0; i < 1500000; i++) {
            p[i] = new Point(ran.nextDouble()*1000, ran.nextDouble()*1000);
        }

        long t1 = System.currentTimeMillis();
        double dis = Point.findNearestPair(p, true);
        long diff = System.currentTimeMillis() - t1;
        System.out.println(dis + " 归并法用时：" + diff + " ms");

        t1 = System.currentTimeMillis();
        dis = Point.findNearestPair(p, false);
        diff = System.currentTimeMillis() - t1;
        System.out.println(dis + " 快排法用时：" + diff + " ms");

        t1 = System.currentTimeMillis();
        dis = Point.findNearestPair(p);
        diff = System.currentTimeMillis() - t1;
        System.out.println(dis + " 划分法用时：" + diff + " ms");

//        使用暴力法验证一下
//        t1 = System.currentTimeMillis();
//        dis = Point.distance(p[0],p[1]);
//        for (int i = 0; i < p.length - 1; i++) {
//            for (int j = i+1; j < p.length; j++) {
//                dis = Math.min(dis, Point.distance(p[j], p[i]));
//            }
//        }
//        diff = System.currentTimeMillis() - t1;
//        System.out.println(dis + "暴力法用时：" + diff + " ms");


    }


}
class Point implements Comparable<Point>{
    private double x;
    private double y;

    Point(double x, double y){
        this.x = x;
        this.y = y;
    }

    public double getX(){
        return x;
    }
    public double getY(){
        return y;
    }

    public void setX(double x) {
         this.x = x;
     }

    public void setY(double y) {
         this.y = y;
     }

    public int compareTo(Point p1){
        if(this.x == p1.x){
            if(this.y == p1.y) return 0;
            else if(this.y > p1.y) return 1;
            else return -1;
        }
        if(this.x > p1.x) return 1;
        return -1;
    }

//     实现hash方法，以便于使用map等数据结构
    public int hashCode(){
        long h, a;
        a = (long) x;
        h = a >>> 32;
        h ^= a;
        a = (long) y;
        h ^= a >>> 32;
        h ^= a;
        return (int) h;
    }

//     静态方法为类方法，静态变量为类变量，一个类的所有对象公用
//     使用时通过类名.类变量名、类名.类方法名调用
//     私有变量在类内可以直接使用对象名.变量名调用，也可以使用类提供的getter方法调用
//     私有方法不能在类外使用，在类外也不能通过对象名.方法名调用。
    public static void show(Point[] points){
        for (int i = 0; i < points.length; i++) {
            System.out.println("( " + points[i].getX() + ", " + points[i].getY() + " )");
        }
    }
    public static void sort(Point[] points){
        Arrays.sort(points);
    }
    public static void sort(Point[] points, int start, int end){
        Arrays.sort(points, start, end);
    }
    public static void sort(Point[] points, Comparator<Point> cmp){
        Arrays.sort(points, cmp);
    }
    public static void sort(Point[] points, int start, int end, Comparator<Point> cmp){
        Arrays.sort(points, start, end, cmp);
    }
    public static Point[] mergeY(Point[] p, Point[] q){
        int length = p.length + q.length;
        Point[] merge = new Point[length];
        int i = 0, j = 0, ind = 0;
        while (i < p.length && j < q.length){
            if(p[i].y < q[j].y){
                merge[ind++] = p[i++];
            }else {
                merge[ind++] = q[j++];
            }
        }
        while (i < p.length){
            merge[ind++] = p[i++];
        }
        while (j < q.length){
            merge[ind++] = q[j++];
        }
        return merge;
    }
    public static double distance(Point p1, Point p2){
        return Math.sqrt(Math.pow(p2.x - p1.x, 2) + Math.pow(p2.y - p1.y, 2));
    }
    public static Point[] filter(Point[] points, int start, int end, double edge, boolean isL){
        int count = 0;
        for (int i = start; i < end; i++) {
            if(isL && points[i].x > edge) count++;
            if(!isL && points[i].x < edge) count++;
        }
        Point[] fil = new Point[count];
        count = 0;
        for (int i = start; i < end; i++) {
            if (isL && points[i].x > edge) fil[count++] = points[i];
            if (!isL && points[i].x < edge) fil[count++] = points[i];
        }
        return fil;
    }
    public static int binSearchY(Point[] points, Point target){
        int index = Arrays.binarySearch(points, target, ((o1, o2) -> (int)(o1.y - o2.y)));
        return index;
    }
//  将sortedByY按照sortedByX划分为两部分
    public static void partition(Point[] sortedByX, Point[] sortedByY, int start, int end){
        //        对x排序数组建立从点到索引的映射
        Map<Point, Integer> index = new HashMap<>();
        for (int i = start; i < end; i++) {
            index.put(sortedByX[i], i);
        }
//        对于按y排序的数组中每一个点，不能像x排序数组那样简单按照mid值分为左右，而是要对其每一个点，
//        都看看这个点如果按照x排序后应该被分在哪里，因此该数组中点的分组要参照points数组（它是按x排的）
        Point[] temp = Arrays.copyOfRange(sortedByY, start, end);
        int mid = start + (end - start) / 2;
        int i = start, j = mid;
        for (int k = 0; k < temp.length; k++) {
            if(index.get(temp[k]) < mid){//如果该点按照x分组时被分到左边，那就放入左边
                sortedByY[i++] = temp[k];
            }else {
                sortedByY[j++] = temp[k];
            }
        }
    }
    public static double findNearestPair(Point[] points){
//        预处理
        Point[] sortY = Arrays.copyOfRange(points, 0, points.length);
        Point.sort(points);
        Point.sort(sortY, new Comparator<Point>() {
            @Override
            public int compare(Point o1, Point o2) {
                if (o1.y > o2.y) return 1;
                if (o1.y < o2.y) return -1;
                return 0;
            }
        });
//        对sortedY划分
        partition(points, sortY, 0, points.length);

        return nearestPair1(points, sortY, 0, points.length);
    }
    public static double findNearestPair(Point[] points, boolean isMergeSort){
//        首先进行排序预处理
        Point.sort(points);
        if(isMergeSort){
            Point[] temp = new Point[points.length];
            System.out.println("mergeSort");;
            return nearestPair(points, temp, 0, points.length);
        }
        System.out.println("quickSort");
        return nearestPair(points, 0, points.length);
    }
    private static double nearestPair(Point[] points, int start, int end){
//        处理[start, end)区间里的点
        int len = end - start;
        int mid = start + len / 2;//偏右
        if(len < 2){
//            少于两个点，就没有计算的必要了
            return Double.MAX_VALUE;
        }
        if(len == 2){
            return Point.distance(points[start], points[start + 1]);
        }
        if(len == 3){
            double d = Point.distance(points[start], points[start + 1]);
            d = Math.min(d, Point.distance(points[start], points[start + 2]));
            d = Math.min(d, Point.distance(points[start + 1], points[start + 2]));
            return d;
        }
//         先求各自的区域中最小的距离
//        System.out.println(start +", "+ mid +", "+ end);
        double minL = nearestPair(points, start, mid);
        double minR = nearestPair(points, mid, end);
        double delta = Math.min(minL, minR);
//        求中线
        double midline = (points[mid - 1].x + points[mid].x) / 2;
//        对于在中线delta范围内的点，再求它们（分别位于两侧）的距离
//         对于左边的点p（x1,y1）,只需要考虑右边范围内的点q（x2,y2），其中y2需要满足：y1-delta < y2 < y1+delta
//         容易证明这样的q至多只有6个

//        分别找出位于两边边界区内的点
        double left = midline - delta, right = midline + delta;
        Point[] P = Point.filter(points, start, mid, left, true);
        Point[] Q = Point.filter(points, mid, end, right, false);
//         对右边界区的点按y坐标排序
        Point.sort(Q, new Comparator<Point>() {
            @Override
            public int compare(Point o1, Point o2) {
                if (o1.y > o2.y) return 1;
                if (o1.y < o2.y) return -1;
                return 0;
            }
        });

        for (int i = 0; i < P.length; i++) {//对于左区域中的每个点
            double down = P[i].y - delta;
            int ind = Point.binSearchY(Q, new Point(0.0, down));
//            ind为负数说明没找到，从-1到-length-1，代表了不同的没找到的位置
//            没找到的位置一共length+1个，正好对应了这些负数（从左往右数）
//            <-1>0<-2>1<-3>2<-4>3<-5>这是二分查找返回的值对应的下标索引
            ind = ind < 0 ? -ind - 1 : ind;
            for (int j = ind; j < 6 + ind && j < Q.length; j++) {//最多只找6次
//                System.out.println("j: " + j + "; length: " + Q.length);
                delta = Math.min(delta, Point.distance(P[i], Q[j]));
            }
        }
        return delta;
    }

    private static double nearestPair(Point[] points, Point[] temp, int start, int end){
//        处理[start, end)区间里的点
        int len = end - start;
        int mid = start + len / 2;//偏右
        if(len < 2){
//            少于两个点，就没有计算的必要了
            temp[start] = points[start];
            return Double.MAX_VALUE;
        }
        if(len == 2){
            if(points[start].y < points[start+1].y){
                temp[start] = points[start];
                temp[start + 1] = points[start + 1];
            }else {
                temp[start] = points[start + 1];
                temp[start + 1] = points[start];
            }
            return Point.distance(points[start], points[start + 1]);
        }
        if(len == 3){
            Point[] p = new Point[1];
            Point[] q = new Point[2];
            p[0] = points[start];
//            q需要有序才能参加归并
            if(points[start + 1].y < points[start+2].y){
                q[0] = points[start + 1];
                q[1] = points[start + 2];
            }else {
                q[1] = points[start + 1];
                q[0] = points[start + 2];
            }

            Point[] t = Point.mergeY(p, q);
            for (int i = 0; i < 3; i++) {
                temp[start + i] = t[i];
            }

            double d = Point.distance(points[start], points[start + 1]);
            d = Math.min(d, Point.distance(points[start], points[start + 2]));
            d = Math.min(d, Point.distance(points[start + 1], points[start + 2]));
            return d;
        }
//        从此开始进入分治递归
//         先求各自的区域中最小的距离
//        System.out.println(start +", "+ mid +", "+ end);
        double minL = nearestPair(points, temp, start, mid);
        double minR = nearestPair(points, temp, mid, end);
        double delta = Math.min(minL, minR);
//        求中线
        double midline = (points[mid - 1].x + points[mid].x) / 2;
//        对于在中线delta范围内的点，再求它们（分别位于两侧）的距离
//         对于左边的点p（x1,y1）,只需要考虑右边范围内的点q（x2,y2），其中y2需要满足：y1-delta < y2 < y1+delta
//         容易证明这样的q至多只有6个

//        分别找出位于两边边界区内的点
        double left = midline - delta, right = midline + delta;
        Point[] P = Point.filter(points, start, mid, left, true);
//        temp中的点是按照y值排序的，其中[start, mid)和[mid, end)中都是分别有序的
        Point[] Q = Point.filter(temp, mid, end, right, false);

        for (int i = 0; i < P.length; i++) {//对于左区域中的每个点
            double down = P[i].y - delta, up = P[i].y + delta;
            int ind = Point.binSearchY(Q, new Point(0.0, down));
//            ind为负数说明没找到，从-1到-length-1，代表了不同的没找到的位置
//            没找到的位置一共length+1个，正好对应了这些负数（从左往右数）
//            <-1>0<-2>1<-3>2<-4>3<-5>这是二分查找返回的值对应的下标索引
            ind = ind < 0 ? -ind - 1 : ind;
            for (int j = ind; j < Q.length && Q[j].y < up; j++) {//最多只找6次
                delta = Math.min(delta, Point.distance(P[i], Q[j]));
            }
        }
//        对temp进行归并排序,需要n次
//        注意这个归并的位置很重要，不能在计算距离之前归并，那样的话，就是左右两个边界区混合在一起了
//        我们只需要右半区有序即可。因为可以看到，在未优化的时候，我们也只是仅对右半区排了序
        Point[] l = new Point[mid - start];
        int idx = 0;
        for (int i = start; i < mid; i++){
            l[idx++] = temp[i];
        }
        Point[] r = new Point[end - mid];
        idx = 0;
        for (int i = mid; i < end; i++) {
            r[idx++] = temp[i];
        }
        Point[] mer = Point.mergeY(l, r);
        for (int i = start; i < end; i++) {
            temp[i] = mer[i - start];
        }
        return delta;
    }

    private static double nearestPair1(Point[] points, Point[] sortedByY, int start, int end){
//        处理[start, end)区间里的点
        int len = end - start;
        int mid = start + len / 2;//偏右
        if(len < 2){
//            少于两个点，就没有计算的必要了
            return Double.MAX_VALUE;
        }
        if(len == 2){
            return Point.distance(points[start], points[start + 1]);
        }
        if(len == 3){
            double d = Point.distance(points[start], points[start + 1]);
            d = Math.min(d, Point.distance(points[start], points[start + 2]));
            d = Math.min(d, Point.distance(points[start + 1], points[start + 2]));
            return d;
        }

//        从此开始进入分治递归
//        进行划分操作
        partition(points, sortedByY, start, end);
//         先求各自的区域中最小的距离
        double minL = nearestPair(points, sortedByY, start, mid);
        double minR = nearestPair(points, sortedByY, mid, end);
        double delta = Math.min(minL, minR);
//        求中线
        double midline = (points[mid - 1].x + points[mid].x) / 2;
//        分别找出位于两边边界区内的点
        double left = midline - delta, right = midline + delta;
        Point[] P = Point.filter(points, start, mid, left, true);
        Point[] Q = Point.filter(sortedByY, mid, end, right, false);

        for (int i = 0; i < P.length; i++) {//对于左区域中的每个点
            double down = P[i].y - delta, up = P[i].y + delta;
            int ind = Point.binSearchY(Q, new Point(0.0, down));
            ind = ind < 0 ? -ind - 1 : ind;
            for (int j = ind; j < Q.length && Q[j].y < up; j++) {//最多只找6次
                delta = Math.min(delta, Point.distance(P[i], Q[j]));
            }
        }
        return delta;
    }
}

6. 结果分析

运行结果：

可见在150万数据时，最快的是快排，归并和划分反而较慢，按照上面的分析好像不太对。我想了一下，可能是因为，快排法不需要额外的空间，但是归并法和划分法都需要额外申请大量空间，在申请空间并赋值的过程消耗了大量时间。