Distance between point and line segments

Question

How would you determine the shortest distance between a point and one or more segments?

For example, what is the shortest distance between the point and the two segments below? Clearly the point is closer to the right line segment. That means the shortest distance between the point and the line segments is equal to the distance between the point and the right segment.

point = {{8, 15}}
lines = {{{20, 10}, {11, 27}}, {{11, 27}, {1, 27}}}
Show[Graphics[{Thick, Line@lines}], 
Graphics[{PointSize[Large], Pink, Point@point}]]

Can you think of a way that also works with an arbitrary number of segments?

Answer 1

In version 10 (now available publicly through the Programming Cloud), use RegionNearest

point = {8, 15};
line = Line[{{{20, 10}, {11, 27}}, {{11, 27}, {1, 27}}}];

np = RegionNearest[line, point]
(* {5663/370, 6981/370} *)

Graphics[{PointSize[Large], {line, Point[point]}, {Red, Point[np]}}]

Answer 2

Using the parametric equation of the lines.

point = {8, 15}
lines = {{{20, 10}, {11, 27}}, {{11, 27}, {1, 27}}}

distance[{start_, end_}, pt_]:= Module[{param = ((pt - start).(end - start))/Norm[end - start]^2},
                               EuclideanDistance[pt, start + Clip[param , {0, 1}] (end - start)]];

Min[distance[#, point] & /@ lines]

(*159/Sqrt[370]*)

Plotting isodistance lines

l  = Table[{Cos[2 Pi n/5], Sin[2 Pi n/5]}, {n, 1, 6}];
l1 = Partition[Riffle[l, l[[2 ;;]]], 2];
Quiet@ContourPlot[Min[distance[#, {x, y}] & /@ l1], {x, -2, 2}, {y, -2, 2}]

Answer 3

Not as elegant as belisarius' first answer but was a bit faster; now slower than the new one. A bit of explanation: If the point on the line segment $PQ$ closest to $X$ is on $PQ$, then the exterior (supplementary) angles of the triangle $PXQ$ at $P$ and $Q$ will be greater than a right angle and the cosines will be negative. [Edit: To clarify, cosAngles is actually a list of Dot products, which are just positive numbers times the cosines; and we only use the sign.]

lines = RandomReal[{0, 10}, {10, 2, 2}];
point = RandomReal[{0, 10}, {2}];

cosAngles = Dot @@@ Partition[Differences[{#1, #2, #3, #1}], 2, 1] &; 
nearest[{P_, Q_}, X_] := Module[{points},
  If[And @@ Negative[cosAngles[X, P, Q]],
   P + Projection[X - P, Q - P], First@Nearest[{P, Q}, X]]
  ];
findNearest[lines_, point_] := 
  Nearest[nearest[#, point] & /@ lines, point];
distance[lines_, point_] := 
  Norm[First@findNearest[lines, point] - point];

Here's the output:

point
findNearest[lines, point]
distance[lines, point]

{7.48559, 3.5353}

{{8.20107, 4.23693}}

1.0021

Here's a picture of what's going on:

g := Module[{distPts, nearPt},
  distPts = nearest[#, point] & /@ lines;
  nearPt = First@Nearest[distPts, point];
  Graphics[{Line@lines, Blue, Thin, Line[{point, #} & /@ distPts], 
    Point[nearest[#, point] & /@ lines], Thick, Line[{point, nearPt}],
     Red, Point@point}]
  ]

g

Same lines, different point. The nearest point is an end point of a segment.

Block[{point = {8, 9}},
 Print[distance[lines, point]];
 g
 ]

1.46937

Answer 4

The numbers shown are the normal distances. 8.266 is the smaller distance. ref above for the equation of the normals.

x0 = 8; y0 = 15;
point = {{x0, y0}};
(*x1=8;y1=20; x2=10;y2=22; x3=12;y3=28*)
x1 = 20; y1 = 10; x2 = 11; y2 = 27; x3 = 1; y3 = 27;
lines = {{{x1, y1}, {x2, y2}}, {{x2, y2}, {x3, y3}}};
v1 = {y2 - y1, -(x2 - x1)};
r1 = {x1 - x0, y1 - y0};
v2 = {y3 - y2, -(x3 - x2)};
r2 = {x3 - x0, y3 - y0};

d1 = Abs@Dot[v1, r1]/Norm[v1];
d2 = Abs@Dot[v2, r2]/Norm[v2];

Labeled[
 Graphics[{
   {Thick, Line@lines}, {PointSize[Large], Pink, Point@point}
   }, Axes -> True, GridLines -> {Range[5, 20, 1], Range[10, 30, 1]}, 
  GridLinesStyle -> LightGray],
 Grid[{{"distance to first line ", N@d1}, 
   {"distance to second line ", N@d2}}, Alignment -> Left]
 ]

Answer 5

This solution I am about to present is effectively a blend of Michael's and belisarius's approaches. To wit,

PointLineDistance[pt_, {s1_, s2_}] :=
                  With[{tp = s1 - pt}, EuclideanDistance[tp, Projection[tp, s2 - s1]]]

segs = {{{20, 10}, {11, 27}}, {{11, 27}, {1, 27}}};
nf = Nearest[segs -> Automatic, DistanceFunction -> PointLineDistance];

PointLineDistance[{8, 15}, Extract[segs, nf[{8, 15}]]]
   159/Sqrt[370]

Trying to reproduce bel's isodistance plot on my box gives something rather different:

segs = Partition[Table[Through[{Cos, Sin}[2 Pi n/5]], {n, 5}] // N, 2, 1, 1];
nf = Nearest[segs -> Automatic, DistanceFunction -> PointLineDistance];

ContourPlot[PointLineDistance[{x, y}, Extract[segs, nf[{x, y}]]], {x, -2, 2}, {y, -2, 2},
            AspectRatio -> Automatic, ColorFunction -> "ThermometerColors"]

Finally, here's a three dimensional test:

BlockRandom[SeedRandom[42, Method -> "MKL"]; (* for reproducibility *)
            segs = RandomReal[{-1, 1}, {50, 2, 3}]];

nf = Nearest[segs -> Automatic, DistanceFunction -> PointLineDistance];

Graphics3D[{Line[segs],
            {Directive[Red, AbsoluteThickness[4]], Line[Extract[segs, nf[{0, 0, 0}]]]},
            {Directive[Red, AbsolutePointSize[6]], Point[{0, 0, 0}]}}]

Answer 6

p1 = {{8, 15}};
polyline = {{{20, 10}, {11, 27}}, {{11, 27}, {1, 27}}};

Discretizing the segments:

(* decrease dt for better result *)
data = With[{dt = .1},
  Flatten[polyline /. r : {{_, _}, {_, _}} :>
     Table[s[[1]] + t (s[[2]] - s[[1]]), {s, r}, {t, 0., 1, dt}], 1]];

p2 = Nearest[data, p][[1]]
(* {{15.5, 18.5}} *)

EuclideanDistance @@ Join[p1, p2]
(* 8.27647 *)

Graphics[{
  Thick, Line@polyline,
  PointSize[Large], Pink, Point[p1~Join~p2]}]

I think this should work in 3D too and with any number of joined or disjoint segments.

Answer 7

Another approach to this problem is to use the image processing functions and operate directly on the image of the lines. Here we import the two lines from the OPs question, and take the DistanceTransform. Finding the value of the distance transform at the orange point gives the distance (in pixels) of that point from that line.

img = Import["https://i.sstatic.net/Inf9t.png"]; 
point = Round[First@ImageValuePositions[img, RGBColor[1, 0.5, 0.5]]]; 
distImg = DistanceTransform[img];
distImg // ImageAdjust
ImageValue[distImg, point]

91.418

Of course one would need to know the relationship between the scaling of the drawing and the distance covered by one pixel in order to convert back to the original units. Chances are this method will be very fast, especially when there are many lines (because the computation time is independent of the number of lines). It also works for curves of any shape, so it is very general. On the other hand, the answer is only approximate, as it works on the rasterized pixel positions and not on the underlying real-valued positions.