Shortest distance between points algorithm

Question

Given a set of points on a plane, find the shortest line segment formed by any two of these points.

How can I do that? The trivial way is obviously to calculate each distance, but I need another algorithm to compare.

Pettifer

http://blogs.mathworks.com/videos/2008/06/02/matlab-puzzler-finding-the-two-closest-points/ This has a whole big bunch of different solutions. Kinda neat to see how many ways there are to the same answer. – Sirmons

There is relatively simple line sweep algorithm that can find the closest pair of points in O(nlogn) time. This article has a short explanation. – Encephalography

Answer 1

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

The problem can be solved in O(n log n) time using the recursive divide and conquer approach, e.g., as follows:

• Sort points along the x-coordinate
• Split the set of points into two equal-sized subsets by a vertical line x = xmid
• Solve the problem recursively in the left and right subsets. This will give the left-side and right-side minimal distances dLmin and dRmin respectively.
• Find the minimal distance dLRmin among the pair of points in which one point lies on the left of the dividing vertical and the second point lies to the right.
• The final answer is the minimum among dLmin, dRmin, and dLRmin.

Oneeyed

You didn't explain how finding dLRmin be done in O(n) time to make the entire algorithm O(nlogn) time – Lewandowski

Would be great to have a code example. I'm not convinced this will evaluate all points that need to be considered. – Tamalatamale

The Wikipedia link is useful, but I would expect some more detailed, or more intuitive explanation, that would add value to the linked article, rather than just copying the content of the article that is not even complete and does not explain a key observation as mentioned also by @KapilD. – Behold

Answer 2

I can't immediately think of a quicker alternative than the brute force technique (although there must be plenty) but whatever algorithm you choose don't calculate the distance between each point. If you need to compare distances just compare the squares of the distances to avoid the expensive and entirely redundant square root.

Thordia

Answer 3

One possibility would be to sort the points by their X coordinates (or the Y -- doesn't really matter which, just be consistent). You can then use that to eliminate comparisons to many of the other points. When you're looking at the distance between point[i] and point[j], if the X distance alone is greater than your current shortest distance, then point[j+1]...point[N] can be eliminated as well (assuming i<j -- if j<i, then it's point[0]...point[i] that are eliminated).

If your points start out as polar coordinates, you can use a variation of the same thing -- sort by distance from the origin, and if the difference in distance from the origin is greater than your current shortest distance, you can eliminate that point, and all the others that are farther from (or closer to) the origin than the one you're currently considering.

Chiffonier

Answer 4

You can extract the closest pair in linear time from the Delaunay triangulation and conversly from Voronoi diagram.

Iodism

This seems promising. Would be nice to have a code example – Tamalatamale

Answer 5

There is a standard algorithm for this problem, here you can find it: http://www.cs.mcgill.ca/~cs251/ClosestPair/ClosestPairPS.html

And here is my implementation of this algo, sorry it's without comments:

    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;
}

Irritating

why is convexHull listed but not used? ccw is also nowhere used. Why the value 100? – Tamalatamale

Answer 6

From your question it is not clear if you are looking for the distance of the segment, or the segment itself. Assuming you are looking for the distance (the segment in then a simple modification, once you know which are the two points whose distance is minimal), given 5 points, numbered from 1 to 5, you need to

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.

If I am not wrong, given the commutativity of the distance you do not need to perform other comparisons. In python, may sound like something

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

That can be tested with something like the following code:

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))

If more than two points can have the same minimal distance, and you are looking for the segments, you need again to modify the proposed code, and the output will be the list of points whose distance is minimal (or couple of points). Hope it helps!

Club

Answer 7

Here is a code example demonstrating how to implement the divide and conquer algorithm. For the algorithm to work, the points x-values must be unique. The non-obvious part of the algorithm is that you must sort both along the x and the y-axis. Otherwise you can't find minimum distances over the split seam in linear time.

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)

Trebuchet