算法导论-第33章-计算几何学（最近点对问题）

最近点对问题

问题描述：在 $n\ge 2$ 个点的集合 $Q$ 中寻找最近点对的问题，“最近”指的是欧几里得距离最小，即点 $p_1=(x_1,y_1)$ 和 $p_2=(x_2,y_2)$ 之间的欧几里得距离 $d=\sqrt{(x_1-x_2{)}^2+(y_1-y_2{)}^2}$ 最小。集合 $Q$ 中的两个点可能会重合，这种情况下，它们之间的距离为0。

现实应用：最近点对问题可以应用于交通控制系统中。为检测出潜在的碰撞事故，在空中或海洋交通控制系统中，需要识别出两个距离最近的交通工具。

---

求平面上$n$个顶点的最近点对问题，对文件data.txt的所有点，求距离最近的点对及其距离。

输入：字符串文件data.txt。每行三个数，分别表示顶点编号和其横纵坐标。

输出：输出距离最近的一对顶点编号，及其距离值。

暴力法

遍历所有点，计算两点之间的距离。如果当前点对之间的距离更小，则更新距离最近的点对信息。

分治法

算法思想：

分治算法的每一次递归调用的输入为子集 $P\subseteq Q$ 以及数组 $X$ 和 $Y$，数组 $X$ 和 $Y$ 均包含 $P$ 中的所有点。按照 $x$ 坐标单调递增的方式对数组 $X$ 中的点排序，按照 $y$ 轴坐标单调递增的方式对数组 $Y$ 中的点排序。

输入为 $P$、$X$ 和 $Y$ 的递归调用首先检查是否有 $|P|\le 3$ 成立。如果成立，则执行暴力法，求 $\left(\begin{array}{c}n\\ 2\end{array}\right)$ 个点对中的最近点对；如果不成立，则递归调用如下分治法。

递归法三步骤：

分解（Divide）：找出一条垂直线 $l$，垂直线 $l$ 满足能够将点集 $P$ 分为两部分 $P_L$ 和 $P_R$ ：使得 $|P_L|=\lceil |P|/2\rceil ,|P_R|=\lfloor |P|/2\rfloor$，即两部分点的个数之差不超过 1。$P_L$ 中所有的点都在直线 $l$ 上或在 $l$ 的左侧，$P_R$ 中所有的点都在直线 $l$ 上或在 $l$ 的右侧。数组 $X$ 被划分为两个数组 $X_L$ 和 $X_R$，分别包含 $P_L$ 和 $P_R$ 中的点，并按 $x$ 轴坐标单调递增的方式排序。类似地，将数组 $Y$ 划分为两个数组 $Y_L$ 和 $Y_R$，分别包含 $P_L$ 和 $P_R$ 中的点，并按 $y$ 轴坐标单调递增的方式排序。

解决（Conquer）：把 $P$ 划分为 $P_L$ 和 $P_R$ 后，再进行两次递归调用，一次找出 $P_L$ 中的最近点对，另一次找出 $P_R$ 中的最近点对。第一次调用的输入为子集 $P_L$、数组 $X_L$ 和 $Y_L$；第二次调用的输入为子集 $P_R$、$X_R$ 和 $Y_R$。令 $P_L$ 和 $P_R$ 返回的最近点对的距离分别为 ${\delta }_L$ 和 ${\delta }_R$，并且置 $\delta =min({\delta }_L,{\delta }_R)$。

合并（Combine）：最近点对要么是某次递归找出距离为 $\delta$ 的点对，要么是 $P_L$ 中的一个点和 $P_R$ 中的一个点组成的点对。如果存在这样的一个点，则点对中的两个点于直线 $l$ 的距离必定都在 $\delta$ 之内。因此，它们必定都处于以直线 $l$ 为中心、宽度为 $2\delta$ 的垂直带型区域内。为了找出这样的点对（如果存在），算法要做如下工作：

1. 建立一个数组 $Y^{\prime }$，它是把数组 $Y$ 中所有不在宽度为 $2\delta$ 的垂直带型区域内的点去掉后所得的数组。数组 $Y^{\prime }$ 与 $Y$ 一样，是按 $y$ 轴坐标排序的。
2. 对数组 $Y^{\prime }$ 中的每一个点 $p$，算法试图找出距离 $p$ 在 $\delta$ 以内的点。在 $Y^{\prime }$中仅需考虑紧随 $p$ 后的7个点。算法计算出从 $p$ 到这7个点的距离，并记录下 $Y^{\prime }$ 的所有点对中最近点对的距离 ${\delta }^{\prime }$。
3. 返回 $\delta$ 和 ${\delta }^{\prime }$ 中较小的那个距离以及对应的点对。

---

代码实现：

代码中封装了两个类，一个是Point，记录的是点的信息，包括编号、$x$ 轴坐标、$y$ 轴坐标；另一个是PointResult，记录的是最近点对信息，包括两个Point和它们之间的距离diastance。

在Point类中，还实现了按照 $x$ 轴坐标和 $y$ 轴坐标对Point对象进行定制排序。

在主类FindTheClosestPointPair中，实现了三个方法，分别是计算欧几里得距离、暴力法求最近点对和分治法求最近点对。

• 计算欧几里得距离：直接按照$d=\sqrt{(x_1-x_2{)}^2+(y_1-y_2{)}^2}$定义计算。
• 暴力法：遍历所有情况，求 $\left(\begin{array}{c}n\\ 2\end{array}\right)$ 个点中的最近点对。
• 分治法：先按 $x$ 轴坐标进行定制排序。点集中点的个数小于或等于3，则直接使用暴力法求解最近点对，保存到PointResult，这也是递归的出口。否则，从中间将点集分为两部分，递归地求解左右两部分。但是要考虑最近点对中的点是一个在垂直线左侧，一个在垂直线上或垂直线右侧的情况。
  • 遍历范围内的点，找出在以垂直线为中心，宽度为2*PointResult.distance范围内的所有点；
  • 按 $y$ 轴坐标进行定制排序；
  • 对于垂直带型区域内的某个点，计算所有和它 $y$ 轴坐标相差不超过PointResult.distance的点距，如果有比之前PointResult中记录小的点距，就更新PointResult。

Java代码实现：

import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Paths;
import java.util.*;

/**
 * 将点的信息（编号、横坐标、纵坐标）封装为Point类
 */
class Point {
    private int number;
    private double x;
    private double y;

    public Point() {
    }

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

    public int getNumber() {
        return number;
    }

    public void setNumber(int number) {
        this.number = number;
    }

    public double getX() {
        return x;
    }

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

    public double getY() {
        return y;
    }

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

    /**
     * 按照x轴坐标单调递增进行定制排序
     * @param points
     * @return Point[]
     */
    public static Point[] compareX(Point[] points) {
        Arrays.sort(points, (o1, o2) -> (int) (o1.getX() - o2.getX()));
        return points;
    }

    /**
     * 按照y轴坐标单调递增定制排序
     * @param points
     * @return Point[]
     */
    public static Point[] compareY(Point[] points) {
        Arrays.sort(points, (o1, o2) -> (int) (o1.getY() - o2.getY()));
        return points;
    }

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

/**
 * 将要返回的结果封装为一个点对类，包括两个点，以及距离
 */
class PointResult{
    private Point p1;
    private Point p2;
    private double distance;

    public PointResult() {
    }

    public PointResult(Point p1, Point p2, double distance) {
        this.p1 = p1;
        this.p2 = p2;
        this.distance = distance;
    }

    public Point getP1() {
        return p1;
    }

    public void setP1(Point p1) {
        this.p1 = p1;
    }

    public Point getP2() {
        return p2;
    }

    public void setP2(Point p2) {
        this.p2 = p2;
    }

    public double getDistance() {
        return distance;
    }

    public void setDistance(double distance) {
        this.distance = distance;
    }

    @Override
    public String toString() {
        return "PointResult{" +
                "p1=" + p1 +
                ", p2=" + p2 +
                ", distance=" + distance +
                '}';
    }
}

public class FindTheClosestPointPair {

    /**
     * 计算两个点之间的欧几里得距离
     * @param p1
     * @param p2
     * @return distance
     */
    public static double distance(Point p1, Point p2) {
        return Math.sqrt((p1.getX() - p2.getX()) * (p1.getX() - p2.getX()) + (p1.getY() - p2.getY()) * (p1.getY() - p2.getY()));
    }

    /**
     * 寻找最近点对（暴力解法）
     * @param points
     * @param left
     * @param right
     * @return PointResult
     */
    public static PointResult violentSolution(Point[] points, int left, int right) {
        PointResult pointResult = new PointResult();
        pointResult.setDistance(Double.MAX_VALUE);
        for (int i = left; i < right; i++) {
            for (int j = i + 1; j < right; j++) {
                double dis = distance(points[i], points[j]);
                if (dis < pointResult.getDistance()) {
                    pointResult.setDistance(dis);
                    pointResult.setP1(points[i]);
                    pointResult.setP2(points[j]);
                }
            }
        }
        return pointResult;
    }

    /**
     * 寻找最近点对（分治法）
     * @param points
     * @param left
     * @param right
     * @return PointResult
     */
    public static PointResult divideAndConquer(Point[] points, int left, int right) {
        PointResult pointResult = null;
        if (right - left <= 3) { // 当点数小于3的时候，直接使用暴力解法求最近点对
            return violentSolution(points, left, right);
        }

        int mid = (left + right) / 2; // 分治（从中点分开）
        PointResult leftPointResult = divideAndConquer(points, left, mid); // 递归计算左部分最近点对
        PointResult rightPointResult = divideAndConquer(points, mid, right); // 递归计算右部分最近点对
        pointResult = leftPointResult.getDistance()  < rightPointResult.getDistance() ? leftPointResult : rightPointResult; // 根据距离判断返回的最近点对是哪个

        // 根据上面求得的distance，找出处于以直线l为中心、宽度为2*distance的垂直带型区域内的点集
        List<Point> area = new ArrayList<>();
        for (int i = left; i < right; i++) {
            if (Math.abs(points[i].getX() - points[mid].getX()) < pointResult.getDistance()) {
                area.add(points[i]);
            }
        }
        // 对上面垂直带型区域内的点集按照y轴坐标定制排序
        Point.compareY(points);

        // 遍历垂直带型区域，要求y轴坐标只差不超过pointResult.distance（因为超过了也就没有必要计算欧几里得距离了，必定比之前得到的resultPoint.distance大）
        for (int i = 0; i < area.size(); i++) {
            for (int j = i + 1; j < area.size() && ((area.get(j).getY() - area.get(i).getY()) < pointResult.getDistance()); j++) {
                double dis = distance(area.get(i), area.get(j));
                if (dis < pointResult.getDistance()) { // 如果有这样的点对，更新pointResult
                    pointResult.setDistance(dis);
                    pointResult.setP1(area.get(i));
                    pointResult.setP2(area.get(j));
                }
            }
        }
        return pointResult;
    }


    public static void main(String[] args) throws IOException {
        // 读取data.txt文件，每一行包含顶点的信息，按空格分割
        String path = "C:\\Projects\\IDEAProjects\\algorithms\\src\\main\\java\\ch33\\data.txt";
        List<String> readAllLines = Files.readAllLines(Paths.get(path));

        Point[] points = new Point[readAllLines.size()];

        for (int i = 0; i < points.length; i++) {
            String[] line = readAllLines.get(i).split("\\s+");
            points[i] = new Point();
            points[i].setNumber(Integer.parseInt(line[0]));
            points[i].setX(Double.parseDouble(line[1]));
            points[i].setY(Double.parseDouble(line[2]));
            //System.out.println(points[i]);
        }

        // 计算算法耗时
        long start = System.currentTimeMillis();

         暴力解法
        //PointResult pointResult = violentSolution(points, 0, points.length);

        // 分治法
        Point.compareX(points);
        PointResult pointResult = divideAndConquer(points, 0, points.length);

        long end = System.currentTimeMillis();

        System.out.println(pointResult);
        System.out.println("算法耗时: " + (end - start) + " ms");
    }
}

结果与分析

在使用朴素算法（暴力）解决最近点对问题的过程中，只需简单的查看所有 $$\begin{gathered} \left(\begin{array}{c}n\\ 2\end{array}\right) \\ =\Theta (n^2) \end{gathered}$$ 个点对。而分治算法，其运行时间可以用 $T(n)=2T(n/2)+\mathrm{O}(n)$ 来描述。因此，该算法的运行时间仅为 $\mathrm{O}(n\log n)$。