点间最短距离算法

Question

给定平面上的一组点，找到由这些点中的任意两个点形成的最短线段。

我该怎么做？最简单的方法显然是计算每个距离，但我需要另一种算法来比较。

Answer 1

http://en.wikipedia.org/wiki/Closest_pair_of_points

这个问题可以用 O( n log n) 使用递归分治方法的时间，例如，如下所示：

• 沿 x 坐标对点进行排序
• 通过垂直线将点集分成两个大小相等的子集 x = xmid
• 在左侧递归地解决问题和右子集。这将分别给出左侧和右侧最小距离 dLmin 和 dRmin。
• 求一对点之间的最小距离 dLRmin，其中一个点位于垂直分割线的左侧，第二个点位于垂直线的右侧。
• 最终的答案是 dLmin、dRmin 和 dLRmin 中的最小值。

Answer 2

我无法立即想到比暴力技术更快的替代方法（尽管一定有很多），但无论您选择什么算法都不会计算每个点之间的距离。如果您需要比较距离，只需比较距离的平方即可，以避免昂贵且完全多余的平方根。

Answer 3

一种可能性是按 X 坐标（或 Y 坐标——哪个并不重要，只要保持一致）对点进行排序。然后，您可以使用它来消除与许多其他点的比较。当您查看点[i]和点[j]之间的距离时，如果单独的X距离大于当前的最短距离，则点[j+1]...点[N]可以被消除为好吧（假设i——如果j，那么点[0]...point[i]被消除）。

如果您的点以极坐标开始，您可以使用同一事物的变体 - 按距原点的距离排序，如果距原点的距离差大于当前的最短距离，您可以消除该点，以及所有其他比您当前正在考虑的距离原点更远（或更近）的位置。

Answer 4

您可以从 Delaunay 三角剖分 中提取线性时间内最接近的对，反之亦然从 Voronoi 图。

Answer 5

这个问题有一个标准算法，在这里你可以找到它：
http://www.cs.mcgill.ca/~cs251/ClosestPair/ ClosestPairPS.html

这是我对该算法的实现，抱歉没有注释：

    static long distSq(Point a, Point b) {
    return ((long) (a.x - b.x) * (long) (a.x - b.x) + (long) (a.y - b.y) * (long) (a.y - b.y));
}

static long ccw(Point p1, Point p2, Point p3) {
    return (long) (p2.x - p1.x) * (long) (p3.y - p1.y) - (long) (p2.y - p1.y) * (long) (p3.x - p1.x);
}

static List<Point> convexHull(List<Point> P) {

    if (P.size() < 3) {
        //WTF
        return null;
    }

    int k = 0;

    for (int i = 0; i < P.size(); i++) {
        if (P.get(i).y < P.get(k).y || (P.get(i).y == P.get(k).y && P.get(i).x < P.get(k).x)) {
            k = i;
        }
    }

    Collections.swap(P, k, P.size() - 1);

    final Point o = P.get(P.size() - 1);
    P.remove(P.size() - 1);


    Collections.sort(P, new Comparator() {

        public int compare(Object o1, Object o2) {
            Point a = (Point) o1;
            Point b = (Point) o2;

            long t1 = (long) (a.y - o.y) * (long) (b.x - o.x) - (long) (a.x - o.x) * (long) (b.y - o.y);

            if (t1 == 0) {
                long tt = distSq(o, a);
                tt -= distSq(o, b);
                if (tt > 0) {
                    return 1;
                } else if (tt < 0) {
                    return -1;
                }
                return 0;

            }
            if (t1 < 0) {
                return -1;
            }
            return 1;

        }
    });



    List<Point> hull = new ArrayList<Point>();
    hull.add(o);
    hull.add(P.get(0));


    for (int i = 1; i < P.size(); i++) {
        while (hull.size() >= 2 &&
                ccw(hull.get(hull.size() - 2), hull.get(hull.size() - 1), P.get(i)) <= 0) {
            hull.remove(hull.size() - 1);
        }
        hull.add(P.get(i));
    }

    return hull;
}

static long nearestPoints(List<Point> P, int l, int r) {


    if (r - l == P.size()) {

        Collections.sort(P, new Comparator() {

            public int compare(Object o1, Object o2) {
                int t = ((Point) o1).x - ((Point) o2).x;
                if (t == 0) {
                    return ((Point) o1).y - ((Point) o2).y;
                }
                return t;
            }
        });
    }

    if (r - l <= 100) {
        long ret = distSq(P.get(l), P.get(l + 1));
        for (int i = l; i < r; i++) {
            for (int j = i + 1; j < r; j++) {
                ret = Math.min(ret, distSq(P.get(i), P.get(j)));
            }
        }
        return ret;

    }

    int c = (l + r) / 2;
    long lD = nearestPoints(P, l, c);
    long lR = nearestPoints(P, c + 1, r);
    long ret = Math.min(lD, lR);
    Set<Point> set = new TreeSet<Point>(new Comparator<Point>() {

        public int compare(Point o1, Point o2) {
            int t = o1.y - o2.y;
            if (t == 0) {
                return o1.x - o2.x;
            }
            return t;
        }
    });
    for (int i = l; i < r; i++) {
        set.add(P.get(i));
    }

    int x = P.get(c).x;

    double theta = Math.sqrt(ret);

    Point[] Q = set.toArray(new Point[0]);
    Point[] T = new Point[Q.length];
    int pos = 0;
    for (int i = 0; i < Q.length; i++) {
        if (Q[i].x - x + 1 > theta) {
            continue;
        }
        T[pos++] = Q[i];
    }

    for (int i = 0; i < pos; i++) {
        for (int j = 1; j < 7 && i + j < pos; j++) {
            ret = Math.min(ret, distSq(T[i], T[j + i]));
        }
    }
    return ret;
}

Answer 6

从您的问题来看，不清楚您是在寻找线段的距离，还是线段本身。假设你正在寻找距离（然后简单修改中的段，一旦你知道哪两个点的距离最小），给定5个点，编号从1到5，你需要

compare 1 with 2,3,4,5, then 
compare 2, with 3,4,5, then 
compare 3 with 4,5, then 
compare 4 with 5.

如果我没记错的话，给定距离的交换律不需要执行其他比较。
在Python中，可能听起来像是

import numpy as np
def find_min_distance_of_a_cloud(cloud):
        """
        Given a cloud of points in the n-dim space, provides the minimal distance.
        :param cloud: list of nX1-d vectors, as ndarray.
        :return:
        """
        dist_min = None
        for i, p_i in enumerate(cloud[:-1]):
            new_dist_min = np.min([np.linalg.norm(p_i - p_j) for p_j in cloud[(i + 1):]])
            if dist_min is None or dist_min > new_dist_min:
                dist_min = new_dist_min

        return dist_min

可以用下面的代码进行测试：

from nose.tools import assert_equal

def test_find_min_distance_of_a_cloud_1pt():
    cloud = [np.array((1, 1, 1)), np.array((0, 0, 0))]
    min_out = find_min_distance_of_a_cloud(cloud)
    assert_equal(min_out, np.sqrt(3))


def test_find_min_distance_of_a_cloud_5pt():
    cloud = [np.array((0, 0, 0)),
             np.array((1, 1, 0)),
             np.array((2, 1, 4)),
             np.array((3, 4, 4)),
             np.array((5, 3, 4))]
    min_out = find_min_distance_of_a_cloud(cloud)
    assert_equal(min_out, np.sqrt(2))

如果两个以上的点可以具有相同的最小距离，并且您正在寻找线段，则需要再次修改建议的代码，并且输出将是距离最小的点（或几个点）的列表。希望有帮助！

Answer 7

这是一个代码示例，演示如何实现分而治之算法。为了使算法正常工作，点的 x 值必须是唯一的。该算法的不明显部分是您必须沿 x 轴和 y 轴排序。否则，您无法在线性时间内找到裂缝上的最小距离。

from collections import namedtuple
from itertools import combinations
from math import sqrt

IxPoint = namedtuple('IxPoint', ['x', 'y', 'i'])
ClosestPair = namedtuple('ClosestPair', ['distance', 'i', 'j'])

def check_distance(cp, p1, p2):
    xd = p1.x - p2.x
    yd = p1.y - p2.y
    dist = sqrt(xd * xd + yd * yd)
    if dist < cp.distance:
        return ClosestPair(dist, p1.i, p2.i)
    return cp

def closest_helper(cp, xs, ys):
    n = len(xs)
    if n <= 3:
        for p1, p2 in combinations(xs, 2):
            cp = check_distance(cp, p1, p2)
        return cp

    # Divide
    mid = n // 2
    mid_x = xs[mid].x
    xs_left = xs[:mid]
    xs_right = xs[mid:]
    ys_left = [p for p in ys if p.x < mid_x]
    ys_right = [p for p in ys if p.x >= mid_x]

    # Conquer
    cp_left = closest_helper(cp, xs_left, ys_left)
    cp_right = closest_helper(cp, xs_right, ys_right)
    if cp_left.distance < cp_right.distance:
        cp = cp_left
    else:
        cp = cp_right

    ys_strip = [p for p in ys if abs(p.x - mid_x) < cp.distance]
    n_strip = len(ys_strip)
    for i in range(n_strip):
        for j in range(i + 1, n_strip):
            p1, p2 = ys_strip[j], ys_strip[i]
            if not p1.y - p2.y < cp.distance:
                break
            cp = check_distance(cp, p1, p2)
    return cp

def closest_pair(points):
    points = [IxPoint(p[0], p[1], i)
              for (i, p) in enumerate(points)]
    xs = sorted(points, key = lambda p: p.x)
    xs = [IxPoint(p.x + i * 1e-8, p.y, p.i)
          for (i, p) in enumerate(xs)]
    ys = sorted(xs, key = lambda p: p.y)
    cp = ClosestPair(float('inf'), -1, -1)
    return closest_helper(cp, xs, ys)