[SOLVED] "Inside polygon" detection and the Illumination Problem

[MarkMLl]

This is not in any way urgent or important, but I thought it might be interesting to throw it open to discussion.

I live a few hundred yards from an ex-WWII airfield now used for private flying, and I rather enjoy the view and the chatter over the RT. However, particularly late in the flying season, some pilots tend to get very sloppy with their approach pattern and overfly areas they should not. The principal offender appears to be a flying instructor, who gives every sign of not imparting the importance of proper procedures to his students.

We're in a shadow as far as ADS-B feeds to e.g. FlightRadar24 are concerned: nothing lower than roughly 1000' (300ish meters) shows up. So over the last few weeks I've been using one of my SDRs with dump1090-mutability (spitting out 100s of Mb per day), and writing Pascal code to filter the output: ignore everything above a certain height, ignore everything outside a defined polygon, generate a ticket file in .gpx format for what's left. There's a lot of subtlety has to go into that, for example many planes have an uncalibrated altimeter feed to their transponder.

I found this paper http://web.archive.org/web/20080812141848/http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ for polygon detection, and have written this code where the input is lat,long expressed as degrees with decimal fraction:

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

1. type
2.   vertex= record
3.             lat, long: double
4.           end;
5.  
6. var
7.   vertices: array of vertex;
8.  
9. (* Working from the current aircraft location, find the angle subtended by each
10.   edge of the polygon by looking at successive pairs of vertices. Assume that
11.   the current location is inside the polygon if the overall angle is 360 degrees.
12. *)
13. function inPoly(here: vertex): boolean;
14.  
15. var
16.   i: integer;
17.   deg1, deg2, subtends: double;
18.   total: double= 0.0;
19.  
20. begin
21.  
22. (* Relative to the fixed point, determine the angle subtended by each edge of   *)
23. (* the polygon. Save this in degrees for relatively easy debugging.             *)
24.  
25.   for i := 0 to Length(vertices) - 2 do begin
26.     deg1 := ArcTan2(vertices[i].lat - here.lat, vertices[i].long - here.long) * (180 / Pi);
27.     deg2 := ArcTan2(vertices[i + 1].lat - here.lat, vertices[i + 1].long - here.long) * (180 / Pi);
28.     subtends := Abs(deg1 - deg2);
29.     if subtends > 180.0 then
30.       subtends := 360.0 - subtends;
31.     total += subtends;
32. // WriteLn(deg1:12:5, deg2:12:5, subtends:12:5);
33.   end;
34. // WriteLn(total:12:5);
35.  
36. (* If the total angle is 360 degrees, then we must be inside the polygon. Any   *)
37. (* smaller angle is outside, but I don't know quite how much slop I should      *)
38. (* allow... testing suggests that the transition is abrupt even without taking  *)
39. (* the sign of the angles into account.                                         *)
40.  
41.   result := (total >= 359.99) and (total < 360.01)
42. // WriteLn('Score: ', total:12:5)
43. end { inPoly } ;
44.  

It actually works surprisingly well: no points result in an angle substantially more than 360 degrees, and there's an extremely rapid falloff below that point.

However, I'm left with a nagging feeling it could be improved by retaining the sign of the angles: which would result in a total angle of 360 degrees if the point were inside the polygon, and zero outside.

The problem is that to get that working, it is obvious that the points that make up the boundary of the polygon would have to be sorted such that progressing round it in a clockwise direction would always result in a +ve angle (or vice-versa, the important thing being consistency).

That is simple enough if the polygon is non-reentrant, at which point the ordering relative to any point inside the boundary would be the same as the ordering relative to the centroid.

However, I have a very unpleasant feeling that for a reentrant polygon this becomes the "Illumination Problem" described by Penrose et al. (https://en.wikipedia.org/wiki/Illumination_problem) and that even if the approach worked it would require sorting the points for every run: possibly relying on a-priori knowledge of whether the test point was inside the polygon.

Are my misgivings well-founded, or is there an easy solution which will allow this algorithm to benefit from retaining the sign?

MarkMLl

[Dzandaa]

Hi MarkMLI

Funny, last week, I wad also playing with dump1090-mutability and Lazarus to display planes near me on a map

For other algo's look at thin link.

https://www.eecs.umich.edu/courses/eecs380/HANDOUTS/PROJ2/InsidePoly.html

B->

[MarkMLl]

Yes, I saw that link and considered some of the other algorithms, but the "angle" one appeared to be a natural solution since I have the vertices rather than (some definition of) the edges.

I'll probably put my stuff up on Github sooner rather than later, but it's looking fairly robust: .gpx output works well with e.g. QGIS.

Of course, on the last decent flying day everybody appeared to be very much on their best behaviour, almost as though Something Had Been Said by one of our less-tolerant neighbours. And that /particularly/ applied to the Persistent Offender...

Mark

ps https://forum.planefinder.net/threads/ads-b-diy-antenna.23/ might interest, but I was getting acceptable results even with a slightly-too-long whip.

[Niglo]

Very simple, very fast, no trigonometry
I've been using it since 30 years!

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

1.   function PointInPolygon(x, y : Single) : Boolean;
2.  
3.  // The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
4.    // with some minor modifications for speed.  It returns 1 for strictly
5.    // interior points, 0 for strictly exterior, and 0 or 1 for points on
6.    // the boundary.
7.   var
8.     I, J: Integer;
9.   begin
10.     Result:=False;                                      
11.      for I:=0 to pred(High(Poly)) do begin
12.       j := succ(i);
13.       if ((((Poly[I].Y<=Y) and (Y<Poly[J].Y)) or
14.         ((Poly[J].Y<=Y) and (Y<Poly[I].Y))) and
15.         (X<(Poly[J].X-Poly[I].X)*(Y-Poly[I].Y)/(Poly[J].Y-Poly[I].Y)+Poly[I].X)) then
16.         Result:=not Result;
17.     end;
18.  end;
19.  

[MarkMLl]

Let me think about that one...

However I must emphasise that I was asking about the peculiarities of the specific algorithm I was using, not for an alternative.

MarkMLl

[MarkMLl]

Using standard optimisation level (so as to be "debugger friendly" etc.)  Wm. Randolph Franklin' algorithm is faster by a very small number of percent (definitely no more than 5%); increasing the optimisation level (without recompiling the RTL) improves things a little. I feel that the performance gain is offset by the loss in clarity.

The program is fairly I/O intensive, but the above stands even if only the "user" time is used. I suspect that the situation would be very different if using the sort of workstations that were available to him when he developed this algorithm, or with an in-memory operation assisted by a GPU or the sort of vectorisation implemented by a Cray etc.

Updated with test results:

Code: Text  [Select][+][-]

1. #!/bin/bash
2.  
3. cat dump1090_1697961664.txt dump1090_1697961664.txt > 2.txt
4. cat 2.txt 2.txt > 4.txt
5. cat 4.txt 4.txt > 8.txt
6. cat 8.txt 8.txt > 16.txt
7.  
8. time ./poly1090-x86_64-linux -- \
9.  50.8982,0.1345 50.9044,0.1522 50.8874,0.1629 50.8810,0.1458 \
10.  16.txt > /dev/null
11.  
12. rm 2.txt 4.txt 8.txt 16.txt
13.  

Original algorithm:

real    0m25.155s
user    0m21.842s
sys     0m1.793s

"Non-trig" algorithm:

real    0m25.923s
user    0m22.861s
sys     0m1.942s

Both above compiled at O4. So there really is nothing in it when working on a fairly large file, which I put down to good integration of the transcendentals with modern x86 implementations.

In a similar vein, doing the calculation with radians rather than degrees /should/ result in an improvement, but in practice it's lost to the I/O overhead.

MarkMLl

[alpine]

*snip*

However, I have a very unpleasant feeling that for a reentrant polygon this becomes the "Illumination Problem" described by Penrose et al. (https://en.wikipedia.org/wiki/Illumination_problem) and that even if the approach worked it would require sorting the points for every run: possibly relying on a-priori knowledge of whether the test point was inside the polygon.

Are my misgivings well-founded, or is there an easy solution which will allow this algorithm to benefit from retaining the sign?

What can be the relation with the "Illumination Problem" actually? I fail to see it. Can you explain a bit more?

Your version won't work with a concave polygon. Compare original version vs yours (by changing {$if 0}):

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

1. uses
2.   math;
3.  
4. type
5.   vertex= record
6.             lat, long: double
7.           end;
8.  
9. var
10.   vertices: array of vertex;
11.  
12. (* Working from the current aircraft location, find the angle subtended by each
13.   edge of the polygon by looking at successive pairs of vertices. Assume that
14.   the current location is inside the polygon if the overall angle is 360 degrees.
15. *)
16. function inPoly(here: vertex): boolean;
17.  
18. var
19.   i: integer;
20.   deg1, deg2, subtends: double;
21.   total: double= 0.0;
22.  
23. begin
24.  
25. (* Relative to the fixed point, determine the angle subtended by each edge of   *)
26. (* the polygon. Save this in degrees for relatively easy debugging.             *)
27.  
28.   for i := 0 to Length(vertices) - 2 do begin
29.     deg1 := ArcTan2(vertices[i].lat - here.lat, vertices[i].long - here.long) * (180 / Pi);
30.     deg2 := ArcTan2(vertices[i + 1].lat - here.lat, vertices[i + 1].long - here.long) * (180 / Pi);
31.  
32.     {$if 0}
33.     subtends := Abs(deg1 - deg2);
34.     if subtends > 180.0 then
35.       subtends := 360.0 - subtends;
36.     {$else}
37.     subtends := deg1 - deg2;
38.     while subtends > 180 do
39.       subtends := subtends - 360;
40.     while subtends < -180 do
41.       subtends := subtends + 360;
42.     {$endif}
43.  
44.     total += subtends;
45. // WriteLn(deg1:12:5, deg2:12:5, subtends:12:5);
46.   end;
47. // WriteLn(total:12:5);
48.  
49. (* If the total angle is 360 degrees, then we must be inside the polygon. Any   *)
50. (* smaller angle is outside, but I don't know quite how much slop I should      *)
51. (* allow... testing suggests that the transition is abrupt even without taking  *)
52. (* the sign of the angles into account.                                         *)
53.  
54.   result := (total >= 359.99) and (total < 360.01)
55. // WriteLn('Score: ', total:12:5)
56. end { inPoly } ;
57.  
58. const
59.   here: vertex = (lat: 25; long: 50 );
60.   r: boolean = false;
61.  
62. begin
63.   SetLength(vertices, 6);
64.   vertices[0].long := 0; vertices[0].lat := 0;
65.   vertices[1].long := 25; vertices[1].lat := 50;
66.   //vertices[2].long := 50; vertices[2].lat := 50;
67.   vertices[2].long := 80; vertices[2].lat := 25;
68.   vertices[3].long := 100; vertices[3].lat := 150;
69.   vertices[4].long := 150; vertices[4].lat := 0;
70.   vertices[5].long := 0; vertices[5].lat := 0;
71.  
72.   r := inPoly(here);
73.  
74. end.

[MarkMLl]

What can be the relation with the "Illumination Problem" actually? I fail to see it. Can you explain a bit more?

The illumination problem, by definition, has a reentrant outline. The result of this is that- I think- if you wanted to order the vertices such that the angles were always +ve as you moved around the edges you'd have to take your position into account.

Thanks for the fix, which I'll check and incorporate. There's a few more cases I want to look at, e.g. a plane moving along the Meridian which is inconveniently close to us.

I put the code on Github earlier since I've got other things looming which might need a clear desk and fast footwork.

https://github.com/MarkMLl/utils1090

MarkMLl

[Mr.Madguy]

Standard way to solve this problem, that is used in games: to turn pairs of points into lines and then use ax + by >= c or ax + by <= c test. It's enough for convex polygon. If polygon isn't convex - BSP-tree should be used instead. May be this can be useful for you.

[alpine]

What can be the relation with the "Illumination Problem" actually? I fail to see it. Can you explain a bit more?

The illumination problem, by definition, has a reentrant outline. The result of this is that- I think- if you wanted to order the vertices such that the angles were always +ve as you moved around the edges you'd have to take your position into account.

Set (unordered) of vertices can give you multiple concave polygons. A classical optimization problem.

[Mr.Madguy]

Well, if polygon isn't convex, then it's vertices have to be sorted, otherwise it would most likely be ambiguous. Convex polygon can be restored though. Via simple rule: for every pair of vertices, that form edge, all other vertices should lie at the same side of this edge.

[MarkMLl]

Thanks for the fix, which I'll check and incorporate. There's a few more cases I want to look at, e.g. a plane moving along the Meridian which is inconveniently close to us.

I can confirm that my original code was broken with even simple reentrant edges, but your fix works very nicely with complex shapes. I can also confirm that the code posted by Niglo works equally well, but there's too much I/O overhead to be able to compare speed realistically.

MarkMLl

[jamie]

I posted a demo of locating the center point of a polygon and it works with complex shapes.

Its here on the group somewhere if you are interested or I could repost it.