# Shortest distance from point to line segment in 3D

This function returns the shortest distance from a 3D point P to a line segment defined by two 3D points A and B. It projects P onto the line, then checks if this projected point lies between A and B. If it does, then the distance is calculated using this projected point.

If the projection is outside the line segment, then the distance is calculated to either point A or B, depending on which is closer.

I wanted to get feedback on if this is a clean solution to the problem.

function d = point_seg_dist(P, A, B)
% Returns distance from a 3D point to a line segment
%
% Inputs
%
%   P : vector
%       Position in space
%
%   A : vector
%       Position of line segment endpoint
%
%   B : vector
%       Position of other line segment endpoint
%
% Outputs
%
%   d : double
%       Distance from point to line segment

AP = P - A; % Vector from A to point
AB = B - A; % Vector from A to B

% Project point onto line
P_line = A + dot(AP, AB) / dot(AB, AB) * AB;

if all(A < P_line) && all(P_line < B)
% The point projected onto the line is in between A and B

% Projection of point onto segment is the same
% as projection of point onto line
P_seg = P_line;
else
% The point projected onto the line is outside of A and B

if all(P_line <= A)
% The projected point is closer to A
P_seg = A;
else
% The projected point is closer to B
P_seg = B;
end
end

d = norm(P - P_seg); % Distance to line segment

end


Here is the implementation of the solution provided by Wolfram Alpha, as was suggested in the comments of the OP.

The implemented equation is:

$$d=\frac{|(\mathbf{x}_2-\mathbf{x}_1)\times(\mathbf{x}_2-\mathbf{x}_1)|}{|\mathbf{x}_1-\mathbf{x}_0|}$$

where

$$x_1=(x_1,y_1,z_1)$$ $$x_2=(x_2,y_2,z_2)$$

are two points on the line and $x_0$ is a third point.

function d = point_seg_dist(x0, x1, x2)

% Make sure we are in 3D
if length(x0) == 2
x1(3) = 0;
x2(3) = 0;
x0(3) = 0;
end

% We want column arrays
if size(x0,1) < size(x0,2)
x0 = P';
x1 = x1';
x2 = x2';
end

d = norm(cross(x2-x1 , x1-x0)) / norm(x2-x1);


Apart from a couple of checks on the input, the calculation is on one line and it does not need projections or checks on the projected point's position.