Finding the perpendicular distance between a vertex and a line

[MarkMLl]

I'm using the code below to find the distance between a vertex and one of the edges making up a polyline, hence to determine which edge is nearest etc.

However after peering at the behaviour for a long time I've realised that what appears to be happening is that the edges are extended past their defined endpoints, so when moving the vertex in the vicinity of a pair of lines which intersect at a shallow angle one finds oneself switching from the nearest one to a straight extension of the other.

Code: Pascal  [Select][+][-]

1. type
2.   TVertexMtrs= record
3.                  lat, lon: double
4.                end;
5.  
6.   TPolylineMtrs= array of TVertexMtrs;
7.  
8.   TPolylinesMtrs= array of TPolylineMtrs;
9.  
10. (* Return the distance, in metres, between the vertex (point) and the indicated
11.   segment of the polyline, where the first segment is indexed zero.
12. *)
13. function ToPolylineMtrs(const v: TVertexMtrs; const p: TPolylineMtrs;
14.                                         segment: integer= 0): TMeters;
15.  
16. (* See https://uknowledge.uky.edu/cgi/viewcontent.cgi?article=1040&context=bae_facpub *)
17. (* hence https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points *)
18. (* of which the fragment below is supposed to be a straight transcript.         *)
19.  
20. begin
21.   Assert(Length(p) >= 2, 'Internal error: in ' + {$I %CURRENTROUTINE%} + '(), polyline does not comprise at least two vertices.');
22.  
23. (* A polyline with two points can only have segment zero and so on.             *)
24.  
25.   Assert(Length(p) - 2 >= segment, 'Internal error : in ' + {$I %CURRENTROUTINE%} + '(), segment outside polyline.');
26.   result := Abs(((p[segment + 1].Lat - p[segment].Lat) * v.Lon) -
27.                 ((p[segment + 1].Lon - p[segment].Lon) * v.Lat) +
28.                 (p[segment + 1].Lon * p[segment].Lat) -
29.                 (p[segment + 1].Lat * p[segment].Lon)) /
30.         Sqrt(Sqr(p[segment + 1].Lat - p[segment].Lat) +
31.                 Sqr(p[segment + 1].Lon - p[segment].Lon))
32. end { ToPolylineMtrs } ;
33.  
34.  

Is there a succinct way of getting the distance to the nearest line, without lines being extended beyond their endpoints?

Any thoughts would be appreciated. My apologies for not posting a compilable example, but this is a facet of a much larger graphical application.

For context, see image at https://forum.lazarus.freepascal.org/index.php/topic,67542.msg523485.html#msg523485 noting that there is a slight kink in the approach at the NE of the airfield. If I traverse a vertex across that kinked pair of legs, the code above behaves erratically as it passes over first one of the legs, then the linear extension of the other.

MarkMLl

[Seenkao]

Найди ближайшие течки каждой линии. Те линии, что будут содержать ближайшую точку, с большой вероятностью будут ближайшими. Ну и перпендикуляры придётся проверить, если линии длинные. Так же можно проверить угол, от проверяемой координаты, если он острый, то вторую проверку делать не надо, если тупой, то проверить надо (на перпендикуляр).

Google translate:
Find the nearest points of each line. Those lines that contain the nearest point are most likely to be the nearest. Well, and the perpendiculars will have to be checked if the lines are long. You can also check the angle from the coordinate being checked, if it is acute, then you do not need to do a second check, if it is obtuse, then you need to check (for the perpendicular).

[MarkMLl]

Nearest vertex would help, but the problem is working out the second one so that the appropriate edge can be selected.

Let's simplify it and stipulate that the edges will be relatively in-line: any deviation will be less than 45 degrees.

However the relative lengths of the edges might differ substantially (the runway will almost certainly be much shorter than the approach or departure leg).

MarkMLl

[Seenkao]

У вас рёбра не конечны? Вы их задаёте формулами, а не координатами? Просто я привык, что работая с графикой, обычно задаются конечные точки.

Есть вариант проверять окружностями. Но тут тоже желательно найти пересечение точки и линии. Вы берёте это пересечение как радиус. И всё что попадает в этот радиус или пересекается с ним, может быть ближе найденной точки. Так мы можем отбросить многие варианты.

Google translate:
Are your edges not finite? Do you define them with formulas, not coordinates? I'm just used to the fact that when working with graphics, end points are usually defined.

There is an option to check with circles. But here, too, it is desirable to find the intersection of a point and a line. You take this intersection as a radius. And everything that falls within this radius or intersects with it can be closer than the found point. This way we can eliminate many options.

[MathMan]

Hi MarkMLI,

Let me see if I got this right

* you want the distance, but only if the base of perpendicular is on the segment of the polyline?

If that is correct then the only solution I know would be to define a line as

1.  L := P1 + r*(P2-P1) // where P1, P2 are the start/end-point of the segment
2. calculate the base of perpendicular on L wrt P3 // call it P4
3. translate P4 to the line format of 1. // t.i. calculate the fitting r
4. the distance (= arg( P4-P3 )) is only valid of 0<=r<=1 in 3. holds

Cheers,
MathMan

[lar]

How about this to determine if the perpendicular between point P and a line segment between points P1 and P2 lies on the segment or not:

1. Let D be the perpendicular distance found with the posted formula, D1 the distance between P and P1, and D2 the distance between P and P2.
2. Compute L1 = sqrt(D1*D1 - D*D) and L2 = sqrt(D2*D2 - D*D).
3. Let L be the distance between P1 and P2.
4. if max(L1, L2) > L then the perpendicular point lies outside the line segment.

[MarkMLl]

Are your edges not finite? Do you define them with formulas, not coordinates? I'm just used to the fact that when working with graphics, end points are usually defined.

Edges are finite, defined by points. The problem appears to be that the expression I've got looks for the perpendicular distance to a line passing through the two vertices, without capping the line at the vertices.

What I could do is extend the edges sideways to make regions of some adequate width, and then (use existing code to) determine which region I'm in before separately looking at the distance.

MarkMLl

[MarkMLl]

If that is correct then the only solution I know would be to define a line as

How about this to determine if the perpendicular between point P and a line segment between points P1 and P2 lies on the segment or not:

Thanks, let me think about those.

MarkMLl

[Thaddy]

NOT mine! but it might give you some ideas:

Code: Pascal  [Select][+][-]

1. type
2.   TPoint = record
3.     X, Y: Double;
4.   end;
5.  
6. function DistanceToSegment(const P1, P2, P3: TPoint; out Distance: Double): Boolean;
7. var
8.   SegmentVec, P3Vec: TPoint;
9.   r, SegmentLengthSq: Double;
10.   Projection: TPoint;
11. begin
12.   // 1. Calculate vector components
13.   SegmentVec.X := P2.X - P1.X;
14.   SegmentVec.Y := P2.Y - P1.Y;
15.  
16.   // Handle zero-length segment edge case
17.   SegmentLengthSq := SegmentVec.X * SegmentVec.X + SegmentVec.Y * SegmentVec.Y;
18.   if SegmentLengthSq < 1e-9 then  // Threshold for floating point precision
19.   begin
20.     Distance := Sqrt(Sqr(P3.X - P1.X) + Sqr(P3.Y - P1.Y));
21.     Result := False;  // Consider segment as point, no valid projection
22.     Exit;
23.   end;
24.  
25.   // 2. Calculate parameter r (projection factor)
26.   P3Vec.X := P3.X - P1.X;
27.   P3Vec.Y := P3.Y - P1.Y;
28.   r := (P3Vec.X * SegmentVec.X + P3Vec.Y * SegmentVec.Y) / SegmentLengthSq;
29.  
30.   // 3. Clamp r to [0,1] range for segment validation
31.   if r < 0.0 then
32.     r := 0.0
33.   else if r > 1.0 then
34.     r := 1.0;
35.  
36.   // 4. Calculate projection point
37.   Projection.X := P1.X + r * SegmentVec.X;
38.   Projection.Y := P1.Y + r * SegmentVec.Y;
39.  
40.   // 5. Calculate distance and validate projection
41.   Distance := Sqrt(Sqr(P3.X - Projection.X) + Sqr(P3.Y - Projection.Y));
42.   Result := (r >= 0.0) and (r <= 1.0);  // True only if projection is ON segment
43. end;

Code: Pascal  [Select][+][-]

1. var
2.   P1, P2, P3: TPoint;
3.   Valid: Boolean;
4.   Dist: Double;
5. begin
6.   // Define points
7.   P1 := TPoint.Create(0, 0);
8.   P2 := TPoint.Create(4, 0);
9.   P3 := TPoint.Create(2, 3);
10.  
11.   // Calculate distance
12.   Valid := DistanceToSegment(P1, P2, P3, Dist);
13.  
14.   if Valid then
15.     WriteLn('Valid distance: ', Dist:0:2)
16.   else
17.     WriteLn('Closest point is segment endpoint. Distance: ', Dist:0:2);
18. end.

Returns True only if the perpendicular base lies between P1 and P2 (0 <= r <= 1)

[munair]

Is this close to what you're looking for?

Code: Pascal  [Select][+][-]

1. type
2.   TVertexMtrs = record
3.     lat, lon: double;
4.   end;
5.  
6.   TPolylineMtrs = array of TVertexMtrs;
7.  
8.   TMeters = double;
9.  
10. function ToPolylineMtrs(const v: TVertexMtrs; const p: TPolylineMtrs; segment: integer = 0): TMeters;
11. var
12.   x, y, x1, y1, x2, y2: double;
13.   dx, dy, t, projX, projY: double;
14.   lenSq: double;
15. begin
16.   Assert(Length(p) >= 2, 'Internal error: Polyline does not comprise at least two vertices.');
17.   Assert(Length(p) - 2 >= segment, 'Internal error: Segment outside polyline.');
18.  
19.   // Segment endpoints (p[segment], p[segment + 1]) and point v
20.   x := v.lon;  y := v.lat;
21.   x1 := p[segment].lon;  y1 := p[segment].lat;
22.   x2 := p[segment + 1].lon;  y2 := p[segment + 1].lat;
23.  
24.   // Vector from start to end of segment
25.   dx := x2 - x1;
26.   dy := y2 - y1;
27.  
28.   // Length squared of the segment
29.   lenSq := Sqr(dx) + Sqr(dy);
30.  
31.   // If segment has zero length, return distance to one of the endpoints
32.   if lenSq = 0 then
33.     Exit(Sqrt(Sqr(x - x1) + Sqr(y - y1)));
34.  
35.   // Compute the projection parameter t
36.   t := ((x - x1) * dx + (y - y1) * dy) / lenSq;
37.  
38.   // Clamp t to [0, 1] to stay within the segment
39.   if t < 0 then
40.   begin
41.     // Closest point is the start of the segment
42.     Exit(Sqrt(Sqr(x - x1) + Sqr(y - y1)));
43.   end
44.   else if t > 1 then
45.   begin
46.     // Closest point is the end of the segment
47.     Exit(Sqrt(Sqr(x - x2) + Sqr(y - y2)));
48.   end;
49.  
50.   // Closest point is within the segment, compute the projected point
51.   projX := x1 + t * dx;
52.   projY := y1 + t * dy;
53.  
54.   // Return distance to the projected point
55.   Result := Sqrt(Sqr(x - projX) + Sqr(y - projY));
56. end;
57.  

Example:

Code: Pascal  [Select][+][-]

1. function FindNearestSegment(const v: TVertexMtrs; const p: TPolylineMtrs; out nearestSegment: integer): TMeters;
2. var
3.   i: integer;
4.   dist, minDist: TMeters;
5. begin
6.   Assert(Length(p) >= 2, 'Polyline too short.');
7.   minDist := Infinity;
8.   nearestSegment := 0;
9.  
10.   for i := 0 to Length(p) - 2 do
11.   begin
12.     dist := ToPolylineMtrs(v, p, i);
13.     if dist < minDist then
14.     begin
15.       minDist := dist;
16.       nearestSegment := i;
17.     end;
18.   end;
19.  
20.   Result := minDist;
21. end;

[MarkMLl]

Is this close to what you're looking for?

Oh absolutely /lovely/: it dropped straight in. I'd been scratching my head over various issues (including aircraft appearing to move to the wrong side of their approach or departure track as shown at the bottom of the window) for several days, and only very gradually realised that I'd misinterpreted the description of the expression I'd "borrowed".

Many, /many/ thanks for that.

MarkMLl