Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

Example: point and two lines

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?

share|improve this question
    
    
@YvesKlett Your link refers to the theory shortest distance. I would like to know how you can calculate distance in practice (with Mathematica of course). –  sjdh Feb 24 '13 at 9:02
    
When you say "arbitrary number of lines", do you mean line, line segment, or ray? –  Silvia Feb 24 '13 at 12:10
    
@Silvia I mean the union of the line segments –  sjdh Feb 24 '13 at 12:35

7 Answers 7

up vote 4 down vote accepted

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

enter image description here

share|improve this answer
2  
OR since the damn Programming cloud has not worked (Server error) since it was released, just use the same RegionNearest in v9 as I did here –  RunnyKine Jun 24 at 15:11

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

Mathematica graphics

share|improve this answer

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

Nearest points on line segments

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

Nearest points on line segments II

share|improve this answer
    
I think my new solution is faster :) –  belisarius Feb 24 '13 at 23:06
1  
@belisarius Touché! Nice work. About twice as fast, or a little better, on a 1000 lines. –  Michael E2 Feb 24 '13 at 23:16
    
Actually 60% faster (1.3/0.8). I really liked the other approach (it seemed more "elegant"), but it was really slow –  belisarius Feb 24 '13 at 23:18
    
@belisarius Well I can't find a different idea that's faster. As a math teacher, I admire your geometric distance function more than using MinValue. It's faster still if you skip Module and insert the param formula directly. –  Michael E2 Feb 25 '13 at 11:28

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

Mathematica graphics

share|improve this answer
    
This works for the line segments given by the OP, but doesn't seem to work in some other case, such as when the line segments are changed to x1 = 8; y1 = 20; x2 = 10; y2 = 22; x3 = 12; y3 = 28; –  Mike Z. Feb 24 '13 at 14:32

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"]

isodistance plot

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

nearest 3D segment test

share|improve this answer
    
I think that nf in the first block finds the segment whose extended line is closest to a point. But the problem is about distance to a segment only, not to its extended line, isn't it? –  BoLe Apr 19 '13 at 9:01
    
How then do you define "distance to a segment"? –  J. M. Apr 19 '13 at 9:02
    
Something like Min[EuclideanDistance[p1 + t (p2 - p1), p3]] with 0 < t < 1. –  BoLe Apr 19 '13 at 9:32
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]}]

line

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

share|improve this answer

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["http://i.stack.imgur.com/Inf9t.png"]; 
point = Round[First@ImageValuePositions[img, RGBColor[1, 0.5, 0.5]]]; 
distImg = DistanceTransform[img];
distImg // ImageAdjust
ImageValue[distImg, point]

enter image description here

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.