分治：平面中的最近点对问题的python实现

最新推荐文章于 2023-12-26 12:05:08 发布

原创最新推荐文章于 2023-12-26 12:05:08 发布 · 5.1k 阅读

48 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#算法 #python #数据结构

分治：平面中的最近点对问题的python实现

前言

分治是算法设计中一种强有力的设计思想，其主要步骤包括了：

将待解决问题递归地划分为若干的子实例
确定最小子问题的解（边界条件）
将子问题的解组合起来，得到原问题的解

本文基于对分治算法基本思想的运用，分析有关平面最近点对问题的分值算法的实现，并且用Python给出具体实现。

问题描述

在平面上存在n个不相重合的点（x_i, y_i）组成的点集，i = 1, 2, … , n, 将任意两个点组成的集合称为一个点对，这两个点之间的距离为这个点对的距离，求这个点集中距离最小的点对。

在这个问题一些其它描述中，允许有重合点的存在，但分析方法大同小异。

问题分析

首先，一种直观的想法是将所有点对的距离求出来取最小值，显然这种算法需要先遍历一遍每个点，对于每个点求出其与其它n - 1个点的距离，时间复杂度为O（n²）。下面我们“严格”地遵循分治算法的基本设计步骤来说明如何将复杂度降为O（nlogn）.

子问题划分&子问题组合

要求解一个点集的最近点对，我们可以利用一条切割线将原点集分为左右两部分，求出左部分的最近点对，再求出右部分的最近点对，最后同时考虑左右部分中比较靠近切割线的点集，求出“中间部分”的最近点对，取三者的最小值。

边界条件

当集合里只有两个点时，返回这两个点之间的距离以及这两个点的坐标。
当集合里只有一个点时，最小距离不存在，返回无穷大（上界），点的坐标返回空集。
由于我们需要的是最短距离，那么我们在比较过程中就要将不符合问题定义的情况舍弃，显然我们需要将这样的情况设置成“返回一个极大值”才能将其在寻找最小值得比较中舍弃掉。

 # 边界条件
    if n == 2:
        print('min_dis = %f' % distance(points_x[0], points_x[1]))
        return [distance(points_x[0], points_x[1]), [points_x[0], points_x[1]]]
    if n == 1:
        return [999, [points_x[0]], points_x[0] # 上界设置为999

代码实现

测试的数据为一个列表，其元素为元组，每个元组包含两个实数作为点的x， y坐标。如：

points_y = [(1, 1), (2, 3), (3, 4), (4, 5), (5, 6), (6, 4), (7, 3), (8, 1), (5, 1), (10, 8),
            (20, 8), (8, 7), (8, 8), (18, 4), (2, 5), (6, 10)]

为了更好地将点集分为大致大小相等的两部分（事实上，不少分治算法在子问题划分上都可以尝试将问题划分为两个大致相等的部分去解决），我们先将点集里的点按x坐标升序排序。

def find_nearest_points(points, n):

我们事先封装好求两点间的距离的函数：

def distance(a, b):
    return ((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2) ** 0.5

这个函数用于求解点集points中的最近点对，n为点的个数。我们规定函数的返回值为一个列表，包含两个元素，第一个元素为最近点对的距离，第二个元素为最近点对。即

return [distance, [(a, b), (c, d)]]

下面我们将点集对半分：

 #  将点集对半分
	left_points = points_x[0:int(n / 2)]
    right_points = points_x[int(n / 2):n + 1]
    
 # 中线横坐标
    lx = (points_x[int(n / 2)][0] + points_x[int(n / 2) - 1][0]) / 2  
    
 #求解左右子问题的解
    d1 = find_nearest_points(left_points, points_y, int(n / 2))
    d2 = find_nearest_points(right_points, points_y, int(n / 2) + (n % 2))

接下来我们需要求解最为棘手的问题，找到两个部分中比较靠近中间部分的最小点对。在刚才的分析中，我们只解决了“只考虑左部分的点”与“只考虑右部分的点”两种子情况，那么如果最近点对中的两个点一个位于左部分，一个位于右部分该怎么办呢？

我们来看下面这幅图：在这里插入图片描述
我们已经能够得知左半部分与右半部分的最近点集以及他们的距离，下面我们设置min_set作为左右两部分中距离较小者的解，需要注意的是min_set同样满足我们规定的解的结构。取min_set中的点对赋值给min_points.

	min_set = min(d1, d2, key=lambda x: x[0])
    min_points = min_set[1]

d = min{左部分最小值，右部分最小值}

在我们已知左右两边最近值为d的情况下，我们只需要像上面这幅图一样，考虑以切割线为中线，边长为2d的矩形区域的点，因为如果原问题的解的点对中，一个点位于左部分，一个点位于右部分，那么这两个点一定在这个区域内，否则它们之间的距离不可能比d还要更小。

在确定了讨论范围后，我们将这个区域的点加入到集合T中。

	T = []
    for i in points_y:
        if lx - min_set[0] <= i[0] <= lx + min_set[0]:
            T.append(i)

上述代码中出现了一个新的点集points_y，这是一个将原点集按纵坐标升序排列的点集，为什么要设置这样一个集合我们马上就会谈到。

我们现在的问题就变成了如何确定集合T中点集中的最近点对，并将它的距离与d作比较。当然我们这里不能求出T中每个点对的距离取最小值，那样的话我们的复杂度又重新变回了O（n²）.为了得出一个在O（n）时间内能都解决问题的复杂度，我们需要找到对于T集合中每个点与其它点比较次数的一个上界。这个上界是否存在呢？我们来看下面这幅图：在这里插入图片描述
在一个d*2d的矩形中，我们考虑一个极端情况，那就是上图的每个红点的位置都放上一个点，这些点之间的距离都不超过d，如果有6个以上的点的话无论如何摆放，都不能满足这个条件。我们可以发现，在这样一个矩形中，最多只能容纳六个相互之间距离不超过d的点（这个结论需要基于点集中不存在相互重合的点），现在我们将这个矩形放在我们的点集T中：在这里插入图片描述
在这样的一个矩形内，我们选取矩形内的每一个点，将选取的点与矩形内的至多5个点进行比较，那么在遍历完所有的点后，如果存在一个距离小于d的点对，那么一定可以找到。那么我们如何选取这5个点呢？一个直观的感受是选取某个点“附近”的五个点，但是我们以怎样的指标选取“附近”呢？如果将与所有点横坐标来确定附近点显然不合理，因为在集合T中我们两个点之间的横坐标已经被限定在2d之内，因此我们可以选取纵坐标最接近选取点的5个点来作比较。

但是这样显然又会出现问题，因为我们如果要找到距某个点纵坐标最近的5个点，那么我们在最差的情况下的时间又会变成O（n²），解决这个问题的方法非常简单，那就是再准备一个以y坐标升序排列的点集points_y，我们将points_y内的点集按顺序加入T的过程中，T也是一个有序点集。

那么我们稍微修改一下函数定义，将points_y包含进来。

def find_nearest_points(points_x, points_y, n):

然后开始执行我们刚才所说的步骤：

	min_dis = min_set[0] #设置min_dis用于储存最小值
    for k in T:
        for i in range(1, 5):
            if T.index(k) + i <= len(T) - 1: # 检测下标是否越界
                if T[T.index(k) + i][1] - k[1] < min_set[0]:
                    dis = distance(k, T[T.index(k) + i])
                    if dis < min_dis: #如果找到了一个比d距离更小的点对，更新距离与点坐标
                        min_dis = dis
                        min_points[0] = k
                        min_points[1] = T[T.index(k) + i]
                else:
                    break
            else:
                break
     return [min(min_dis, min_set[0]), min_points]

在调用这个函数之前，需要准备两个原点集以x升序与以y升序排列的集合。
以下是完整代码：

import copy


# 求a，b两点间的距离
def distance(a, b):
    return ((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2) ** 0.5


def find_nearest_points(points_x, points_y, n):

    # 边界条件
    if n == 2:
        return [distance(points_x[0], points_x[1]), [points_x[0], points_x[1]]]
    if n == 1:
        return [999, [points_x[0]], points_x[0]]
    
    #  将点集对半分
    left_points = points_x[0:int(n / 2)]
    right_points = points_x[int(n / 2):n + 1]
    lx = (points_x[int(n / 2)][0] + points_x[int(n / 2) - 1][0]) / 2  # 中线横坐标

    d1 = find_nearest_points(left_points, points_y, int(n / 2))
    d2 = find_nearest_points(right_points, points_y, int(n / 2) + (n % 2))

    min_set = min(d1, d2, key=lambda x: x[0])
    min_points = min_set[1]
   
    # print("d = %f" % min_set[0])

    # print('lx = %f' % lx)
    T = []

    for i in points_y:
        if lx - min_set[0] <= i[0] <= lx + min_set[0]:
            T.append(i)
    # 计算T1中每个点到最多5个点的距离
    min_dis = min_set[0]
    for k in T:
        for i in range(1, 5):
            if T.index(k) + i <= len(T) - 1:
                if T[T.index(k) + i][1] - k[1] < min_set[0]:
                    dis = distance(k, T[T.index(k) + i])
                    if dis < min_dis:
                        min_dis = dis
                        min_points[0] = k
                        min_points[1] = T[T.index(k) + i]
                else:
                    break
            else:
                break
    return [min(min_dis, min_set[0]), min_points]