Puzzle 23 - Distance to line segment

[gmmastros]

Calculate the shortest distance from a point to a line segment.






Given the coordinates of the line segment (shown in blue), calculate the shortest distance to a given point (also shown in blue). The distance that you are to calculate is represented by the red line. If the point is 'outside' the line segment, meaning, you cannot draw a perpendicular line to it, then you should return the distance from the given point to the closest end point (shown in the second image).

The calculations are relatively simple. No points will be rewarded for incorrect distances (within reasonable rounding issues). Max points will be awarded for correct results and fast performance.

[ChibaPet]

Here's mine. I'm bad at math, so I didn't use any of those scary geometry things, except for a distance formula. I'd love to have implemented this without sqrt or exponentiation, but I'm not sure how.

perlLine number Off | Show | Select all

1. #!/usr/bin/perl
2.  
3. use strict;
4. use warnings;
5.  
6. distance_to_line_segment(100, 100, 200, 200, 285, 19);
7. distance_to_line_segment(100, 100, 200, 200, 195, 255);
8.  
9. sub distance_to_line_segment
10. {
11.     my ($ax, $ay, $bx, $by, $cx, $cy) = @_;
12.     my ($c1x, $c1y, $c2x, $c2y, $c3x, $c3y);
13.  
14.     # Find point c1, which is c shifted such that c2x == ax and c2y mirrors
15.     # c1y, across the line segment ax,ay bx,by.
16.     my $slopey = ($ay - $by) / ($ax - $bx);
17.     my $deltax = $cx - $ax;
18.     my $deltay = $deltax * $slopey;
19.     $c1x = $ax;
20.     $c1y = $cy + $deltax * $slopey;
21.  
22.     # Find point c2, which is c shifted such that c2y == ay and c2x mirrors
23.     # c1x, across the line segment ax,ay bx,by.
24.     my $slopex = 1 / $slopey;
25.     $deltay = $cy - $ay;
26.     $deltax = $deltay * $slopex;
27.     $c2x = $cx + $deltax;
28.     $c2y = $ay;
29.  
30.     # Find the midpoint of the line segment c1x,c1y c2x,c2y.
31.     $c3x = ($c1x + $c2x) / 2;
32.     $c3y = ($c1y + $c2y) / 2;
33.  
34.     my $distance;
35.     if ($c3x < $ax and $c3y < $ay) {
36.         # The midpoint of c1x,c1y c2x,c2y is outside point a.
37.         $distance =  sqrt((($cx - $ax) ** 2) + (($cy - $ay) ** 2));
38.     } elsif ($c3x > $bx and $c3y > $by) {
39.         # The midpoint of c1x,c1y c2x,c2y is outside point b.
40.         $distance =  sqrt((($cx - $bx) ** 2) + (($cy - $by) ** 2));
41.     } else {
42.         # The midpoint of c1x,c1y c2x,c2y is inside points a and b.
43.         $distance =  sqrt((($cx - $c3x) ** 2) + (($cy - $c3y) ** 2));
44.     }
45.  
46.     print "Distance to ($ax,$ay $bx,$by) from $cx,$cy is $distance.\n";
47. }
48.  

Code is hidden, SHOW

[ChibaPet]

Here's mine. I'm bad at math, so I didn't use any of those scary geometry things, except for a distance formula. I'd love to have implemented this without sqrt or exponentiation, but I'm not sure how.

Hm, yeah, my solution fails if the points defining the line segment are listed in the reverse order. I'm not going to bother fixing it, but I should figure out which of a and b is closest to the origin, and have that determine which endpoint is checked against "<" and which is checked against ">". I also feel like there are common cases I'm missing, but maybe that would suffice.

Alright, I'll do it, or it's going to gnaw at me.

perlLine number Off | Show | Select all

1. #!/usr/bin/perl
2.  
3. use strict;
4. use warnings;
5.  
6. distance_to_line_segment(100, 100, 200, 200, 285, 19);
7. distance_to_line_segment(100, 100, 200, 200, 195, 255);
8.  
9. sub distance_to_line_segment
10. {
11.     my ($ax, $ay, $bx, $by, $cx, $cy) = @_;
12.     my ($c1x, $c1y, $c2x, $c2y, $c3x, $c3y);
13.  
14.     # Order our points so that the one closer to the origin is considered first.
15.     if (distance($ax, $ay, 0, 0) > distance($bx, $by, 0, 0)) {
16.         my ($tempx, $tempy);
17.  
18.         $tempx = $ax; $tempy = $ay;
19.         $ax = $bx;    $ay = $by;
20.         $bx = $tempx; $by = $tempy;
21.     }
22.  
23.     # Find point c1, which is c shifted such that c2x == ax and c2y mirrors
24.     # c1y, across the line segment ax,ay bx,by.
25.     my $slopey = ($ay - $by) / ($ax - $bx);
26.     my $deltax = $cx - $ax;
27.     my $deltay = $deltax * $slopey;
28.     $c1x = $ax;
29.     $c1y = $cy + $deltax * $slopey;
30.  
31.     # Find point c2, which is c shifted such that c2y == ay and c2x mirrors
32.     # c1x, across the line segment ax,ay bx,by.
33.     my $slopex = 1 / $slopey;
34.     $deltay = $cy - $ay;
35.     $deltax = $deltay * $slopex;
36.     $c2x = $cx + $deltax;
37.     $c2y = $ay;
38.  
39.     # Find the midpoint of the line segment c1x,c1y c2x,c2y.
40.     $c3x = ($c1x + $c2x) / 2;
41.     $c3y = ($c1y + $c2y) / 2;
42.  
43.     my $distance;
44.     if ($c3x < $ax and $c3y < $ay) {
45.         # The midpoint of c1x,c1y c2x,c2y is outside point a.
46.         $distance =  distance($cx, $cy, $ax, $ay);
47.     } elsif ($c3x > $bx and $c3y > $by) {
48.         # The midpoint of c1x,c1y c2x,c2y is outside point b.
49.         $distance =  distance($cx, $cy, $bx, $by);
50.     } else {
51.         # The midpoint of c1x,c1y c2x,c2y is inside points a and b.
52.         $distance =  distance($cx, $cy, $c3x, $c3y);
53.     }
54.  
55.     print "Distance to ($ax,$ay $bx,$by) from $cx,$cy is $distance.\n";
56. }
57.  
58. sub distance
59. {
60.     my ($ax, $ay, $bx, $by) = @_;
61.  
62.     return sqrt((($ax - $bx) ** 2) + (($ay - $by) ** 2));
63. }
64.  

Code is hidden, SHOW

There. Now the order of points in the line segment doesn't matter, at the expense of two more distance calculations per problem.

[AlexCuse]

Nice work ChibaPet, and Welcome

I can't test it, but it looks about like what I had in mind (which is probably bad news for you!)

[koorogi]

I followed a link here from reddit and figured I'd give it a go. Here's my solution in C.

cLine number Off | Show | Select all

1. /* needed for sqrtf and fabsf */
2. #ifndef _ISOC99_SOURCE
3. #define _ISOC99_SOURCE
4. #endif
5.  
6. #include <stdlib.h>
7. #include <stdio.h>
8. #include <math.h>
9.  
10. #define MIN(a,b)                ((a) < (b) ? (a) : (b))
11. #define MAX(a,b)                ((a) > (b) ? (a) : (b))
12.  
13. #define IN_INTERVAL(lo,hi,val)  ((val) >= MIN(lo,hi) && (val) <= MAX(lo,hi))
14.  
15. #define SQUARE(a)               ((a)*(a))
16. #define DIST(x1,y1,x2,y2)       sqrtf(SQUARE(x2-x1) + SQUARE(y2-y1))
17.  
18. /* return distance to the nearest endpoint */
19. static float endpoints(float x1, float y1, float x2, float y2, float px, float py)
20. {
21.     return MIN(DIST(x1,y1, px,py), DIST(x2,y2, px,py));
22. }
23.  
24. /* return distance assuming a horizontal line */
25. static float hline_distance(float x1, float x2, float y, float px, float py)
26. {
27.     if(!IN_INTERVAL(x1,x2, px))
28.         return endpoints(x1,y,x2,y,px,py);
29.     else
30.         return fabsf(py-y);
31. }
32.  
33. /* returns the distance between (px,py) and the line segment from (x1,y1) to (x2,y2) */
34. static float distance(float x1, float y1, float x2, float y2, float px, float py)
35. {
36.     float slope1, slope2, yoff1, yoff2, x, y;
37.  
38.     /* special case horizontal lines */
39.     if(y1 == y2)
40.         return hline_distance(x1, x2, y1, px, py);
41.     /* and vertical lines */
42.     else if(x1 == x2)
43.         return hline_distance(y1, y2, x1, py, px);
44.  
45.     /* for all other lines:
46.      * find the intersection between the line segment and the perpendicular
47.      * line containing the point (px,py).  Return the distance between this
48.      * intersection point and (px,py) if it really lies in the interval,
49.      * otherwise just return the distance to the closest endpoint.
50.      */
51.     slope1 = (y2-y1) / (x2-x1);
52.     slope2 = -1 / slope1;
53.  
54.     yoff1  = y1 - slope1*x1;
55.     yoff2  = py - slope2*px;
56.  
57.     x = (yoff2 - yoff1) / (slope1 - slope2);
58.  
59.     if(!IN_INTERVAL(x1,x2, x))
60.         return endpoints(x1,y1,x2,y2,px,py);
61.  
62.     y = slope1*x + yoff1;
63.     return DIST(x,y, px,py);
64. }
65.  
66. int main(int argc, char **argv)
67. {
68.     if(argc >= 7)
69.         printf("%f\n", distance(atof(argv[1]), atof(argv[2]),
70.                                 atof(argv[3]), atof(argv[4]),
71.                                 atof(argv[5]), atof(argv[6])));
72.     return 0;
73. }
74.  

Code is hidden, SHOW

[chrissie1]

Here's mine. I'm bad at math, so I didn't use any of those scary geometry things, except for a distance formula. I'd love to have implemented this without sqrt or exponentiation, but I'm not sure how.

Sorry to sound stupid but what language is that you are using.

[damber]

Here's mine. I'm bad at math, so I didn't use any of those scary geometry things, except for a distance formula. I'd love to have implemented this without sqrt or exponentiation, but I'm not sure how.

Sorry to sound stupid but what language is that you are using.

ummm... I think my sarcasm filter is a bit wonky... are you serious? (I really can't tell for some reason )

...or didn't you read the first line?

[chrissie1]

Here's mine. I'm bad at math, so I didn't use any of those scary geometry things, except for a distance formula. I'd love to have implemented this without sqrt or exponentiation, but I'm not sure how.

Sorry to sound stupid but what language is that you are using.

ummm... I think my sarcasm filter is a bit wonky... are you serious? (I really can't tell for some reason )

...or didn't you read the first line?

ah math is the language. Didn't know that was a programming language. Or am I misstaken.

My limited English tell me math is short fo mathematics.

ah now I see it's perl.

Nope I missed that one.

I changed the hidecode tags to reflect that for other people like me.

[gmmastros]

I've created an executable to perform this calculation and also to show the task graphically. This should run on any windows 2000+ operating system.

http://tinyurl.com/b66rfr

[AlexCuse]

Chrissie, hopefully there aren't too many people like you

This was much trickier than I expected, but I learned a lot! Cut out some stuff in the LineSegment class (I hope to eventually make this able to determine the intersection point and create the intersecting line as well - perhaps this could be a puzzle for another week?) hopefully didn't break anything

csharpLine number Off | Show | Select all

1. using System;
2.  
3. namespace PointToLine
4. {
5.     class Program
6.     {
7.         static void Main(string[] args)
8.         {
9.             while (true)
10.             {
11.                 Point[] points = GetPoints();
12.                 LineSegment line = new LineSegment() { P1 = points[0], P2 = points[1] };
13.                 Console.WriteLine("Distance is: {0}", PointToLineDistance(points[2], line));
14.                 Console.WriteLine("\n\nWould you like to try again? (Y/N)");
15.                 if (Console.ReadLine().ToUpper() != "Y")
16.                     break;
17.             }
18.             Console.WriteLine("\nGoodbye!");
19.         }
20.  
21.         static Double PointToLineDistance(Point point, LineSegment target)
22.         {
23.             if ((point - target.P2) * (target.P2 - target.P1) > 0)
24.             {
25.                 return Math.Sqrt((target.P2 - point) * (target.P2 - point));
26.             }
27.             else if ((point - target.P1) * (target.P1 - target.P2) > 0)
28.             {
29.                 return Math.Sqrt((target.P1 - point) * (target.P1 - point));
30.             }
31.             else
32.             {
33.                 return Math.Abs(((target.P2 - target.P1) ^ (point - target.P1)) / Math.Sqrt((target.P2 - target.P1) * (target.P2 - target.P1)));
34.             }
35.         }
36.         static Point[] GetPoints()
37.         {
38.             Point[] output = new Point[3];
39.  
40.             Console.WriteLine("Lets define the points - P1 and P2 will define the line segment, P3 will be the point to get distance from...");
41.             int pX, pY;
42.             for (int i = 0; i < 3; i++)
43.             {
44.                 Console.WriteLine("P{0} x: ", i + 1);
45.                 pX = HandleInput();
46.                 Console.WriteLine("P{0} y: ", i + 1);
47.                 pY = HandleInput();
48.                 output[i] = new Point() { X = pX, Y = pY };
49.             }
50.             return output;
51.         }
52.  
53.         static Int32 HandleInput()
54.         {
55.             Int32 x;
56.             while (!Int32.TryParse(Console.ReadLine(), out x))
57.             {
58.                 Console.WriteLine("C'mon man, give me an integer!");
59.             }
60.             return x;
61.         }
62.     }
63.  
64.     class Point
65.     {
66.         public Double X { get; set; }
67.         public Double Y { get; set; }
68.         public override String ToString()
69.         {
70.             return String.Format("X:{0} Y:{1}", X, Y);
71.         }
72.  
73.         /// <summary>Dot Product</summary>
74.         public static double operator *(Point left, Point right)
75.         {
76.             return (left.X * right.X) + (left.Y * right.Y);
77.         }
78.  
79.         /// <summary>Cross Product</summary>
80.         public static double operator ^(Point left, Point right)
81.         {
82.             return (left.X * right.Y) - (left.Y * right.X);
83.         }
84.  
85.         public static Point operator -(Point left, Point right)
86.         {
87.             return new Point() { X = left.X - right.X, Y = left.Y - right.Y };
88.         }
89.  
90.         public static Point operator +(Point left, Point right)
91.         {
92.             return new Point() { X = left.X + right.X, Y = left.Y + right.Y };
93.         }
94.  
95.     }
96.  
97.     class LineSegment
98.     {
99.         public Point P1 { get; set; }
100.         public Point P2 { get; set; }
101.     }
102. }

Code is hidden, SHOW

[gmmastros]

I'm trying to test the code submitted.

So far, I've only tested koorogi's code. His code produces the correct answer and is only slightly slower than my own code, which I will post later.

When I compile my code, it runs in 160 nano-seconds. Koorogi's code (converted to VB, sorry) runs in 300 nano-seconds.

Please, can everyone test their own performance. Run the code a couple million times and post the execution time (divided by the number of loops of course).

[AlexCuse]

A tick is supposedly one hundred nanoseconds, I'm sure its' a bit more than that though. I'm too lazy to do anything more complicated.

These times are on a pretty slow laptop, with a lot of other stuff running.

[hide]Lets define the points - P1 and P2 will define the line segment, P3 will be the
point to get distance from...
P1 x:
100
P1 y:
100
P2 x:
200
P2 y:
200
P3 x:
350
P3 y:
795
Distance is: 613.616329639295
Duration is: 2590000 Ticks or 1.295 Ticks per execution


Would you like to try again? (Y/N)[/hide]

[gmmastros]

Chibapet,

The code you posted works if the line is not horizontal or vertical.

[ChibaPet]

Chibapet,

The code you posted works if the line is not horizontal or vertical.

Yar. Broken.

The following corrects that. It could search for point C being on the line which contains the line segment, but I guess that would only be worthwhile if it's a common case in the data.

Anyway, it's less pretty now, but more correct. Lucky I don't program spacecraft controls, I guess.

perlLine number Off | Show | Select all

1. #!/usr/bin/perl
2.  
3. use strict;
4. use warnings;
5.  
6. distance_to_line_segment(100, 100, 200, 200, 285, 19);
7. distance_to_line_segment(100, 100, 200, 200, 195, 255);
8.  
9. distance_to_line_segment(10, 10, 10, 20, 10, 30);
10. distance_to_line_segment(10, 10, 20, 10, 30, 10);
11. distance_to_line_segment(20, 10, 30, 10, 10, 10);
12. distance_to_line_segment(20, 10, 30, 10, 25, 15);
13.  
14. sub distance_to_line_segment
15. {
16.     my ($ax, $ay, $bx, $by, $cx, $cy) = @_;
17.     my ($c1x, $c1y, $c2x, $c2y, $c3x, $c3y);
18.     my $distance;
19.  
20.     # Order our points so that the one closer to the origin is considered
21.     # first.
22.     if (distance($ax, $ay, 0, 0) > distance($bx, $by, 0, 0)) {
23.         my ($tempx, $tempy);
24.         $tempx = $ax; $tempy = $ay;
25.         $ax = $bx;    $ay = $by;
26.         $bx = $tempx; $by = $tempy;
27.     }
28.  
29.     if ($ax == $bx) {
30.         # Line segment is vertical.
31.         if ($cy < $ay) {
32.             # c is inside of a
33.             $distance =  distance($cx, $cy, $ax, $ay);
34.         } elsif ($cy > $by) {
35.             # c is outside of b
36.             $distance =  distance($cx, $cy, $bx, $by);
37.         } else {
38.             # c is bracketed by a and b
39.             $distance = abs($cx - $ax);
40.         }
41.     } elsif ($ay == $by) {
42.         # Line segment is horizontal.
43.         if ($cx < $ax) {
44.             # c is inside of a
45.             $distance =  distance($cx, $cy, $ax, $ay);
46.         } elsif ($cx > $bx) {
47.             # c is outside of b
48.             $distance =  distance($cx, $cy, $bx, $by);
49.         } else {
50.             # c is bracketed by a and b
51.             $distance = abs($cy - $ay);
52.         }
53.     } else {
54.         # Find point c1, which is c shifted such that c2x == ax and c2y
55.         # mirrors c1y, across the line segment ax,ay bx,by.
56.         my $slopey = ($ay - $by) / ($ax - $bx);
57.         my $deltax = $cx - $ax;
58.         my $deltay = $deltax * $slopey;
59.         $c1x = $ax;
60.         $c1y = $cy + $deltax * $slopey;
61.  
62.         # Find point c2, which is c shifted such that c2y == ay and c2x
63.         # mirrors c1x, across the line segment ax,ay bx,by.
64.         my $slopex = 1 / $slopey;
65.         $deltay = $cy - $ay;
66.         $deltax = $deltay * $slopex;
67.         $c2x = $cx + $deltax;
68.         $c2y = $ay;
69.  
70.         # Find the midpoint of the line segment c1x,c1y c2x,c2y.
71.         $c3x = ($c1x + $c2x) / 2;
72.         $c3y = ($c1y + $c2y) / 2;
73.  
74.         # What is the relation of c3 to ax,ay bx,by?
75.         if ($c3x < $ax and $c3y < $ay) {
76.             # The midpoint of c1x,c1y c2x,c2y is outside point a.
77.             $distance =  distance($cx, $cy, $ax, $ay);
78.         } elsif ($c3x > $bx and $c3y > $by) {
79.             # The midpoint of c1x,c1y c2x,c2y is outside point b.
80.             $distance =  distance($cx, $cy, $bx, $by);
81.         } else {
82.             # The midpoint of c1x,c1y c2x,c2y is inside points a and b.
83.             $distance =  distance($cx, $cy, $c3x, $c3y);
84.         }
85.     }
86.  
87.     print "Distance to ($ax,$ay $bx,$by) from $cx,$cy is $distance.\n";
88. }
89.  
90. sub distance
91. {
92.     my ($ax, $ay, $bx, $by) = @_;
93.  
94.     return sqrt((($ax - $bx) ** 2) + (($ay - $by) ** 2));
95. }
96.  

Code is hidden, SHOW

[gmmastros]

Here is my algorithm:

VBLine number Off | Show | Select all

1. Public Function PerpendicularDistanceToLine(ByVal xPoint As Single, ByVal yPoint As Single, ByVal xA As Single, ByVal yA As Single, ByVal xB As Single, ByVal yB As Single) As Single
2.    
3.     Dim Slope As Single
4.    
5.     Dim xInt As Single
6.     Dim yInt As Single
7.    
8.     If xA = xB Then
9.        
10.         xInt = xA
11.         yInt = yPoint
12.        
13.     ElseIf yA = yB Then
14.        
15.         xInt = xPoint
16.         yInt = yA
17.    
18.     Else
19.        
20.         Slope = ((yB - yA) / (xB - xA))
21.         xInt = (yPoint - yA + xA * Slope + xPoint / Slope) / (1 / Slope + Slope)
22.         yInt = (xPoint - xA + yA * Slope + yPoint / Slope) / (1 / Slope + Slope)
23.    
24.     End If
25.    
26.     If xA < xB Then
27.         If xInt < xA Then
28.             xInt = xA
29.         ElseIf xInt > xB Then
30.             xInt = xB
31.         End If
32.     ElseIf xA > xB Then
33.         If xInt < xB Then
34.             xInt = xB
35.         ElseIf xInt > xA Then
36.             xInt = xA
37.         End If
38.     End If
39.    
40.     If yA < yB Then
41.         If yInt < yA Then
42.             yInt = yA
43.         ElseIf yInt > yB Then
44.             yInt = yB
45.         End If
46.     ElseIf yB < yA Then
47.         If yInt < yB Then
48.             yInt = yB
49.         ElseIf yInt > yA Then
50.             yInt = yA
51.         End If
52.     End If
53.  
54.     PerpendicularDistanceToLine = Sqr((xInt - xPoint) * (xInt - xPoint) + (yInt - yPoint) * (yInt - yPoint))
55.    
56. End Function
57.  

Code is hidden, SHOW

My code has many similarities to code already posted. The main difference is... I calculate the point where the 'distance' line intersects the line segment. This could be anywhere on the line or at one of the end points. From there... it's a simple distance calculation.

[gmmastros]

By the way.... I'd love to see someone post a javascript version. Anyone?

[ChibaPet]

Here is my algorithm:

Hm, I'm not understanding what you're doing with the slope... Math is not my strong point. But does that work? Comments in the code would be a welcome addition!

Also, near the end, I think you want the square root in your distance calculation... Is Sqr the square root in VB? A quick Google doesn't enlighten me overmuch about a VB built-in square root function, but the principle of least surprise says "Sqr" squares its argument.

[gmmastros]

Is Sqr the square root in VB?.

Yes.

http://msdn.microsoft.com/en-us/library/aa263367(VS.60).aspx

I'll write a mathematical proof soon. So that (hopefully), it will be clear how this works. I wrote this code (actually, a slower, uglier version) several years ago, so I don't remember exactly how it was derived. If I recall correctly, I started with the equation for a couple of the lines, and then reduced the formulas a lot to get them in to this form. I tweaked them a lot to get them as fast as possible.

[gmmastros]

Entering equations is not the 'most fun' thing to do, so I wrote it out long hand. I hope nobody has any problem reading this.

Xa, Ya is one point of the line segment and Xb,Yb is the other point. Xp,Yp is the point to calculate the distance for. Xi,Yi is the point on the line where the perpendicular 'distance' line intersect the line segment. In the drawing, it intersects within the end points, but it may end up being outside the end points. This mathematical proof does not take that in to account, but the VB code does.

Also.... shown below is the proof for the X intersect point. Calculations for the Y intersect point is not shown, but is sufficiently similar that is *probably* doesn't need to be shown.



You'll notice that this matches the equation I use in my algorithm for the x-intercept (where the distance line intersects the line segment). Immediately after the calculations, I check the X and Y values to make sure they are on the line. If not, I set them to the closest point.

Then, at the end, I calculate the distance using the normal distance calculation.

[gmmastros]

For what it's worth, I found a bug in my code. If the line segment is actually a point (distance = 0, x and y's are the same), my code calculates the wrong distance, which should actually be a simple point to point distance. Anyway, here's the corrected code.

vbLine number Off | Hide | Select all

1. Public Function PerpendicularDistanceToLine(ByVal xPoint As Single, ByVal yPoint As Single, ByVal xA As Single, ByVal yA As Single, ByVal xB As Single, ByVal yB As Single) As Single
2.    
3.     Dim Slope As Single
4.    
5.     Dim xInt As Single
6.     Dim yInt As Single
7.    
8.     If xA = xB Then
9.        
10.         xInt = xA
11.         yInt = yPoint
12.        
13.     ElseIf yA = yB Then
14.        
15.         xInt = xPoint
16.         yInt = yA
17.    
18.     Else
19.        
20.         Slope = ((yB - yA) / (xB - xA))
21.         xInt = (yPoint - yA + xA * Slope + xPoint / Slope) / (1 / Slope + Slope)
22.         yInt = (xPoint - xA + yA * Slope + yPoint / Slope) / (1 / Slope + Slope)
23.    
24.     End If
25.    
26.     If xA < xB Then
27.         If xInt < xA Then
28.             xInt = xA
29.         ElseIf xInt > xB Then
30.             xInt = xB
31.         End If
32.     ElseIf xA > xB Then
33.         If xInt < xB Then
34.             xInt = xB
35.         ElseIf xInt > xA Then
36.             xInt = xA
37.         End If
38.     ElseIf xA = xB Then
39.         xInt = xA
40.     End If
41.    
42.     If yA < yB Then
43.         If yInt < yA Then
44.             yInt = yA
45.         ElseIf yInt > yB Then
46.             yInt = yB
47.         End If
48.     ElseIf yB < yA Then
49.         If yInt < yB Then
50.             yInt = yB
51.         ElseIf yInt > yA Then
52.             yInt = yA
53.         End If
54.     ElseIf yA = yB Then
55.         yInt = yA
56.     End If
57.  
58.     PerpendicularDistanceToLine = Sqr((xInt - xPoint) * (xInt - xPoint) + (yInt - yPoint) * (yInt - yPoint))
59.    
60. End Function
61.  

The problem was right at the end where I check to see if the xInt is within xA and xB and is yInt is between yA and Yb. Notice the If/Else If blocks. I needed to add a special case where they are equal so that the intersect points would BOTH get set properly.

[ChibaPet]

For what it's worth, I found a bug in my code. If the line segment is actually a point (distance = 0, x and y's are the same), my code calculates the wrong distance, which should actually be a simple point to point distance.

Just tested - mine seems to get that right. I still need to time it, though.

Okay. It's consistently taking five seconds to run a million times. I'm using the original problem, so it's running close to the longest code path - no horizontal or vertical shortcuts. No point reordering, but that's trivial anyway.

I'll get some microsecond resolution soon, and I'll get an average time per run.

[Emtucifor]

May I suggest that

ElseIf yA = yB Then

can just be

Else

[natinja]

A solution based on a Python transcription of the Point class found in Squeak (http://squeak.org).

pythonLine number Off | Show | Select all

1. #!/usr/bin/python
2.  
3. from math import sqrt
4.  
5. class Point:
6.     def __init__(self, x, y):
7.         self.x = x
8.         self.y = y
9.  
10.     def __repr__(self):
11.         """Textual representation of the receiver"""
12.         return "Point(%d, %d)" % (self.x, self.y)
13.  
14.     def __le__(self, other):
15.         return self.x <= other.x and self.y <= other.y
16.  
17.     def __lt__(self, other):
18.         return self.x < other.x and self.y < other.y
19.  
20.     def adhere(self, p1, p2):
21.         """If the receiver lies outside the rectangle formed by p1 and p2, return the nearest point on the boundary of the rectangle, otherwise return self."""
22.         if (p1 <= self and self < p2):
23.             return self
24.         else:
25.             return Point(min(max(p1.x, self.x), p2.x), min(max(p1.y, self.y), p2.y))
26.  
27.     def distance(self, another_point):
28.         """Answer the distance between another_point and the receiver."""
29.         dx = another_point.x - self.x
30.         dy = another_point.y - self.y
31.         return sqrt((dx * dx) + (dy * dy))
32.  
33.     def nearest_point_along_line(self, p1, p2):
34.         """Find the nearest point along a line going from p1 to p2"""
35.         if (p1.x == p2.x):
36.             return Point(p1.x, self.y)
37.         if (p1.y == p2.y):
38.             return Point(self.x, p1.y)
39.         x1 = float(p1.x)
40.         y1 = float(p1.y)
41.         x2 = float(p2.x) - x1
42.         y2 = float(p2.y) - y1
43.         t = (((self.y - y1) / x2) + ((self.x - x1) / y2)) / ((x2 / y2) + (y2 / x2))
44.         return Point(x1 + (x2 * t), y1 + (y2 * t))
45.  
46.     def nearest_point_on_line(self, p1, p2):
47.         return self.nearest_point_along_line(p1, p2).adhere(p1, p2)
48.  
49. # case 1
50. p1 = Point(100, 100)
51. p2 = Point(200, 200)
52. p3 = Point(285, 19)
53. p4 = p3.nearest_point_on_line(p1, p2)
54. print "%s -> %s: %f" % (p3, p4, p3.distance(p4))
55.  
56. # case 2
57. p5 = Point(195, 255)
58. p6 = p5.nearest_point_on_line(p1, p2)
59. print "%s -> %s: %f" % (p5, p6, p5.distance(p6))

Code is hidden, SHOW

Trying to solve the same problem in Squeak is much easier since it already provides a Point class as part of the standard graphic primitives with all the methods for doing the calculation. Note the nice shorthand to instantiate a Point by sending the #@ message and the Y coordinate to a number.

smalltalkLine number Off | Show | Select all

1. ((p := 285@19) nearestPointOnLineFrom: 100@100 to: 200@200) dist: p

Code is hidden, SHOW

[GothGargoyle]

Here's my solution in Java. It actually contains two solutions - the entire thing can basically be done with a built-in function, so as well as using that I went ahead and did the manual calculations for fun

javaLine number Off | Show | Select all

1. package lessthandot;
2.  
3. import java.awt.geom.Line2D;
4. import java.awt.geom.Point2D;
5.  
6. public class DistanceToLine {
7.  
8.     public static void main(final String[] args) {
9.         if (args.length < 6) {
10.             System.out.println("Input should be coordinates of line AB and point C as follows: Ax Ay Bx By Cx Cy");
11.             return;
12.         }
13.         DistanceToLine distanceToLine = new DistanceToLine();
14.        
15.         // Calculate distance using java.awt.geom functions
16.         double distance = distanceToLine.getDistance(Double.valueOf(args[0]), Double.valueOf(args[1]),
17.                 Double.valueOf(args[2]), Double.valueOf(args[3]), Double.valueOf(args[4]), Double.valueOf(args[5]));
18.         System.out.println("Distance from point to line using java.awt.geom is " + distance);
19.        
20.         // Calculate distance manually
21.         distance = distanceToLine.getDistanceManually(Double.valueOf(args[0]), Double.valueOf(args[1]),
22.                 Double.valueOf(args[2]), Double.valueOf(args[3]), Double.valueOf(args[4]), Double.valueOf(args[5]));
23.         System.out.println("Distance from point to line calculated manually is " + distance);
24.     }
25.    
26.     public double getDistance(final double x1, final double y1, final double x2, final double y2, final double px, final double py) {
27.         Point2D point = new Point2D.Double(px, py);
28.         Line2D line = new Line2D.Double(x1, y1, x2, y2);
29.         return line.ptSegDist(point);
30.     }
31.    
32.     public double getDistanceManually(final double x1, final double y1, final double x2, final double y2, final double px, final double py) {
33.         if ((x2 - x1) * (px - x2) + (y2 - y1) * (py - y2) > 0) {
34.             // Point is closest to point B of line AB
35.             return Math.sqrt((x2 - px) * (x2 - px) + (y2 - py) * (y2 - py));
36.         } else if ((x1 - x2) * (px - x1) + (y1 - y2) * (py - y1) > 0) {
37.             // Point is closest to point A of line AB
38.             return Math.sqrt((x1 - px) * (x1 - px) + (y1 - py) * (y1 - py));
39.         } else {
40.             double cross = (x2 - x1) * (px - x1) - (y2 - y1) * (py - y1);
41.             double lineLength = Math.sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
42.             return Math.abs(cross / lineLength);
43.         }
44.     }
45.    
46. }
47.  

Code is hidden, SHOW

[damber]

it would be interesting to see the performance of the builtin function vs the various manual methods here... if anyone has the time to run them through