Efficient method for finding closest point to a line segment from a set of points including line segment's vertices

Question

Let's say I have a list of points (in my case, point objects in a Python implementation). I then have a line segment connecting two of these points.

I want to know if there's a way to efficiently find the point from the list that is closest to the line segment. I realize I could step through all of the points and check distances individually, then select the smallest distance, but I'm eventually hoping to scale the program up to deal with perhaps millions of points. So, you know, maybe something more efficient than that?

If it helps, in terms of CCW all the points being questioned should either be to the "left"/"below" or the "right"/"above" of the line segment; I don't think my implementation will involve checking points to both sides of the segment.

The points will be point objects with (x,y) coordinates and some other junk not directly relevant to this question. The line segments are currently implemented as objects containing references to two point objects (its vertices) and its length.

As I mentioned this is part of a Python implementation. I'm trying to design a way to find a concave hull to a set of points (given some predefined parameters of how to decide whether a point not on the convex hull is on the concave hull or not). I want to find the convex hull first, then for each line segment on the convex hull find the closest point to it, make a triangle with that point, then decide if that triangle bears "deleting such that the internal point is now on the hull.

I considered putting this in Mathematics, but I don't need a distance equation - I need help with an efficient algorithm for finding points closest to a line segment. Note also that I am not looking for the closest point on a line to an input point; I'm looking for the closest point from a set to an input line segment.

Thank you all!

Answer 1

You're going to have to loop through all the points and calculate the distance. There is a great question on StackOverflow about how to calculate the distance: Shortest distance between a point and a line segment

Some of the work can be precalculated, given that you have to do this more than once for a given line segment. You also don't need to figure out the smallest distance, only the smallest distance squared (since squaring is strictly increasing). I converted the top voted answer in that question to a Java version that does this precalculation (in the constructor), and is easier to read and follow. Unfortunately I don't know Python well enough to give you Python code, but you should be able to figure it out.

import java.util.Arrays;
import java.util.List;

public class LineSegment {
    public static class Point {
        public final double x;
        public final double y;

        public Point(double x, double y) {
            this.x = x;
            this.y = y;
        }

        public String toString() {
            return "(" + x + "," + y + ")";
        }
    }

    public static void main(String[] args) {
        LineSegment segment = new LineSegment(new Point(0, 3), new Point(2, 0));
        List<Point> pointList =
                Arrays.asList(new Point[] { new Point(-5, 3), new Point(1, 1),
                        new Point(2, 3), new Point(0, 5) });

        Point answer = segment.closestPoint(pointList);
        System.out.println("The closest point is: " + answer);
    }

    private static double sqr(double x) {
        return x * x;
    }

    private static double distanceSquared(Point v, Point w) {
        return sqr(v.x - w.x) + sqr(v.y - w.y);
    }

    private final Point firstSegPoint;
    private final Point secondSegPoint;
    private final double segmentDistance;
    private double xDifference;
    private double yDifference;

    public LineSegment(Point firstSegPoint, Point secondSegPoint) {
        this.firstSegPoint = firstSegPoint;
        this.secondSegPoint = secondSegPoint;
        this.segmentDistance = distanceSquared(firstSegPoint, secondSegPoint);
        this.xDifference = secondSegPoint.x - firstSegPoint.x;
        this.yDifference = secondSegPoint.y - firstSegPoint.y;
    }

    public Point closestPoint(List<Point> pointList) {
        double minDistance = Double.POSITIVE_INFINITY;
        Point answer = null;

        for (Point point : pointList) {
            double distSquared = distToSegmentSquared(point);
            if (distSquared < minDistance) {
                answer = point;
                minDistance = distSquared;
            }
        }

        return answer;
    }

    private double distToSegmentSquared(Point input) {
        if (segmentDistance == 0)
            return distanceSquared(input, firstSegPoint);

        double xComponent = (input.x - firstSegPoint.x) * xDifference;
        double yComponent = (input.y - firstSegPoint.y) * yDifference;
        double t = (xComponent + yComponent) / segmentDistance;
        if (closestPointIsFirst(t))
            return distanceSquared(input, firstSegPoint);
        if (closestPointIsSecond(t))
            return distanceSquared(input, secondSegPoint);
        Point closestPointOnLine =
                new Point(firstSegPoint.x + t * xDifference, firstSegPoint.y
                        + t * yDifference);
        return distanceSquared(input, closestPointOnLine);
    }

    private boolean closestPointIsFirst(double t) {
        return t < 0;
    }

    private boolean closestPointIsSecond(double t) {
        return t > 1;
    }
}

See the full implementation here: http://ideone.com/fBFwda

Answer 2

For what it is worth, for a static set of points to check the line segment against:

Precomputation using voronoi diagrams (http://en.wikipedia.org/wiki/Voronoi_diagram) is helpful as it reduces the nearest neighbor problem to a polygon line intersection query.

The voronoi diagram is made up of polygons. Each point from the set of points you are trying to compare against is inside it's own polygon. No point is contained within more than 1 polygon and no other points are inside that polygon.

If you only intersect 1 polygon, then you have your answer. If you intersect more than 1 polygon, you must check each polygon you intersected.

Precomputation of the voronoi diagram will not help for dynamic sets of points to query against.

Answer 3

What a fun question! I love grids personally. I love em, I love em, I love em. I wanna partition everything into a grid for no reason.

But the way I'd give this a try would start with the grid. You can partition the points into the grid cells. Each grid cell can be a 32-bit index. Each point can have a next index pointing to the next point in the cell, like so:

Then you can rasterize your line into the grid using Bresenham, like so:

Now to find the closest point to the segment, look for the points inside the cells plotted by the line. If you don't find any points, then expand and look for all points in the neighboring cells and so forth in a breadth-first fashion.

This is probably not the most optimal algorithm and the worst-case scenario is when you have to do this kind of breadth-first search for many phases before finding a single point, but it's the type I always found very easy to tune and you can play with the appropriate density/coarseness of each grid cell (smaller vs. bigger cells) based on your inputs. If you can afford a postprocess after partitioning the points into the cells, you can do a cache-friendliness pass to make sure each point in each cell is contiguous to the next point in the cell for spatial locality.

I'm assuming besides having a boatload of points, you also have many line segments as well to search against. I also used this type of solution except in 3-dimensions and my goal wasn't to find the nearest point to a segment, but which points were closest to which segments. My goal was to find them for 3d line segments representing a bone in a skeletal hierarchy, like so:

Using a 3D bresenham rasterizer using a voxel kind of approach:

... to then produce bone weights for characters:

... and start animating 'em:

In a stress test I was able to generate the skinning weights for a mesh with a million vertices (points) against a few hundred bones (3d line segments) in about 20ms just using this plain old dumb grid approach to it. I wanted that kind of performance so that artists could just draw the bones into the character and start animating it, tweak the bone positions, etc. without waiting for a progress bar each time.

In competing products I often found this process can take minutes and I'm willing to bet they're using much fancier and smarter algorithms than a plain old grid. I'm actually very incompetent mathematically and algorithmically but have managed to outdo colleagues a million times smarter than me using dumb solutions like this just tuned a lot and playing to the idea of bits and bytes instead of a KD-tree with a cutting-edge median splitting algorithm. Sometimes a dumber solution like this can be just fine for your needs provided it's efficient in terms of memory access patterns and so forth as opposed to much fancier solutions like BVHs, KD-trees, quad-trees, octrees. Of course I don't recommend writing a raytracer using a 3-dimensional grid, there are limits, but you can do a lot of things competitively fast using just a plain old grid.

What sort of application does this have in 2-dimensions with the convex hull and everything? I kinda wanna implement it now and see how well the grid fares in such a case.