193. Shortest distance from a point to a line

A recent forum thread asked how to find the shortest distance from a point to a line. As usual, our resident mathematician was there to provide the answer.

In fact, two answers.

A long one,

-- Returns the distance from the point p to the line segment from a to b
-- p,a,b are all vec2s
function lineDist(p,a,b)
    -- Vector along the direction of the line from a to b
    local v = b - a
    -- Project p onto the vector v and compare with the projections of a and b
    if v:dot(p) < v:dot(a) then     -- The projection of p is lower than a, so return the distance from p to a         return p:dist(a)     elseif v:dot(p) > v:dot(b) then
    -- The projection of p is beyond b, so return the distance from p to b
        return p:dist(b)
    end
    -- The projection of p is between a and b, so we need to take the dot product with the orthogonal vector to v, so we rotate v and normalise it.
    v = v:rotate90():normalize()
    -- This is the distance from v to the line segment.
    return math.abs(v:dot(p)-v:dot(a))
end

and a shorter version (one line of code).

function lineDistB(p,a,b)
    return (p-a):dist(math.max(math.min((p-a):dot(b-a)/(b-a):lenSqr(),1),0)*(b-a))
end

plus a bonus function for getting the intersection point

function lineIntersectPoint(p,a,b)
    return math.max(math.min((p-a):dot(b-a)/(b-a):lenSqr(),1),0)*(b-a)+a
end

Here is his explanation of the solution. Any blue text is my comment.

Dot products are the heart of geometry. For this particular problem, the key piece of information is that if you have a unit vector, say u, then the orthogonal projection of a point, say p, onto the line in the direction of u is (p.u) u. If v is a vector which is not necessarily a unit vector, the formula is (p.v) v / (v.v).

So in your code, you first work out the vector along the line (this is the dx = x2 - x1 dy = y2 - y1 part).

Then you work out the square of the length of this vector (temp = dx * dx + dy * dy); this is the same as the dot product of (dx,dy) with itself.

Next, you work out the projection of (px,py) onto the line in the direction of (dx,dy) but as our line doesn’t pass through the origin, we have to adjust it slightly by projecting (px,py) - (x1,y1) so that everything is relative to one end of the line (this is the a=((px-x1)*dx+(py-y1)*dy)/temp part). So in dot product terms, your variable a contains (p.v) / (v.v).

Since we’re dealing with a line segment, we test whether this a is between 0 and 1 and clamp it to that range (I’d probably do a a = math.max(math.min(a,1),0) instead of the conditionals here).

Once you’ve clamped, you compute the point (x1,y1) + a (dx,dy) which is the closest point to (px,py) on the line segment.

Finally, you return its distance from (px,py).

Mathematically, you compute

||p - Max(Min(p.(b - a)/(b-a).(b-a),1),0)(b-a)||

lf you didn’t understand all that, well, nor did I, but I filed the function away in case I ever need it, and I thought I’d share it, because forum threads get buried very quickly.

 

Advertisement