给定平面上的一组点,找出其中任意两个点形成的最短线段。
我怎样才能做到这一点?简单的方法显然是计算每个距离,但我需要另一种算法来比较。
给定平面上的一组点,找出其中任意两个点形成的最短线段。
我怎样才能做到这一点?简单的方法显然是计算每个距离,但我需要另一种算法来比较。
http://en.wikipedia.org/wiki/Closest_pair_of_points
可以使用递归分治法在 O(n log n) 时间内解决该问题,例如,如下所示:
我无法立即想到比蛮力技术更快的替代方法(尽管必须有很多),但是您选择的任何算法都不会计算每个点之间的距离。如果您需要比较距离,只需比较距离的平方,以避免昂贵且完全冗余的平方根。
一种可能性是通过点的 X 坐标(或 Y ——实际上并不重要,只要保持一致)对点进行排序。然后,您可以使用它来消除与许多其他点的比较。当您查看 point[i] 和 point[j] 之间的距离时,如果仅 X 距离大于您当前的最短距离,则可以将 point[j+1]...point[N] 消除为好吧(假设i<j
-- if j<i
,那么它是 point[0]...point[i] 被消除)。
如果您的点从极坐标开始,您可以使用同一事物的变体 - 按与原点的距离排序,如果与原点的距离差大于您当前的最短距离,您可以消除该点,以及所有其他比您当前考虑的更远(或更接近)原点的其他人。
您可以从Delaunay 三角剖分中提取最接近的线性时间对,反之则从Voronoi 图中提取。
这个问题有一个标准算法,在这里你可以找到它: 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;
}
从您的问题来看,不清楚您是在寻找段的距离还是段本身。假设您正在寻找距离(然后进行简单修改,一旦您知道哪些是距离最小的两个点),给定 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))
如果两个以上的点可以具有相同的最小距离,并且您正在寻找线段,则需要再次修改建议的代码,输出将是距离最小的点列表(或几个点)。希望能帮助到你!
这是一个演示如何实现分治算法的代码示例。为了使算法起作用,点 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)