Алгоритм кратчайшего расстояния между точками

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

Проблема может быть решена за время O (n log n) с использованием рекурсивного подхода «разделяй и властвуй», например, следующим образом:

• Сортировать точки по координате x
• Разделите набор точек на два подмножества равного размера вертикальной линией x = xmid
• Решите задачу рекурсивно в левом и правом подмножествах. Это даст минимальные расстояния dLmin и dRmin для левой и правой стороны соответственно.
• Найдите минимальное расстояние dLRmin среди пары точек, в которой одна точка лежит слева от разделяющей вертикали, а вторая - справа.
• Окончательный ответ - это минимум среди dLmin, dRmin и dLRmin.

Я не могу сразу придумать более быструю альтернативу, чем метод грубой силы (хотя их должно быть много), но какой бы алгоритм вы ни выбрали, не вычисляйте расстояние между каждой точкой. Если вам нужно сравнить расстояния, просто сравните квадраты расстояний, чтобы избежать дорогостоящего и полностью избыточного квадратного корня.

Одна из возможностей - отсортировать точки по их координатам X (или Y - на самом деле не имеет значения, какие, просто будьте последовательны). Затем вы можете использовать это, чтобы исключить сравнения со многими другими точками. Когда вы смотрите на расстояние между точкой [i] и точкой [j], если только расстояние X больше, чем ваше текущее кратчайшее расстояние, тогда точка [j + 1] ... точка [N] может быть удалена как хорошо (при условии i<j - если j<i, то удаляются точки [0] ... точка [i]).

Если ваши точки начинаются как полярные координаты, вы можете использовать один и тот же вариант - отсортировать по расстоянию от начала координат, и если разница в расстоянии от начала координат больше, чем ваше текущее кратчайшее расстояние, вы можете удалить эту точку, и все остальные, которые находятся дальше (или ближе) от источника, чем тот, который вы сейчас рассматриваете.

Вы можете извлечь ближайшую пару за линейное время из триангуляции Делоне и с помощью диаграмма Вороного.

Существует стандартный алгоритм решения этой проблемы, вы можете найти его здесь: 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)