Deprecated: Array and string offset access syntax with curly braces is deprecated in /homepages/21/d38531796/htdocs/jose/smfforum/Sources/Subs.php on line 3825
Print Page - Point in polygon

Theo's Forum

IT-Consultant: Charles Pegge => OxygenBasic => Topic started by: Nicola_Piano on April 07, 2022, 04:53:09 PM

Title: Point in polygon
Post by: Nicola_Piano on April 07, 2022, 04:53:09 PM
Hello,
I am looking for an algorithm that allows me to understand if a point, of which I know the coordinates, is inside or outside a polygon of which I know the coordinates of the vertices.
On the internet I found this ...

The code below is from Wm. Randolph Franklin <w...@ecse.rpi.edu>
with some minor modifications for speed. It returns 1 for strictly
interior points, 0 for strictly exterior, and 0 or 1 for points on the boundary. The boundary behavior is complex but determined;
in particular, for a partition of a region into polygons, each point is "in" exactly one polygon. See the references below for more detail.

Code: [Select]
int pnpoly(int npol, float *xp, float *yp, float x, float y)
{
int i, j, c = 0;
for (i = 0, j = npol-1; i < npol; j = i++) {
if ((((yp[i]<=y) && (y<yp[j])) ||
((yp[j]<=y) && (y<yp[i]))) &&
(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))

c = !c;
}
return c;
}

How to translate it into O2?
Cheers
Title: Re: Point in polygon
Post by: Charles Pegge on April 07, 2022, 08:23:20 PM
Hi Nicola,

Here is a partially readable translation:

Code: [Select]
function pnpoly(int npol, float *xp, float *yp, float x, float y) as int
=======================================================================
indexbase 0
int i, j, c = 0;
for i=1 to i<npol
  j=i-1
  'for (i = 0, j = npol-1; i < npol; j = i++) {
  if ((((yp[i]<=y) && (y<yp[j])) ||
    ((yp[j]<=y) && (y<yp[i]))) &&
    (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
    c = not c;
  endif
next
return c;
end function

'TEST
float px={0,1,1,0}
float py={0,0,1,1}
float x,y
int n=4
float x=0.5, y=0.5
print pnpoly(n,px,py,x,y)
print pnpoly(n,px,py,x+0.6,y)
print pnpoly(n,px,py,x,y+0.6)
Title: Re: Point in polygon
Post by: Nicola_Piano on April 07, 2022, 10:13:15 PM
Thanks Charles,
tomorrow I'll try it and tell you.
Title: Re: Point in polygon
Post by: Charles Pegge on April 08, 2022, 01:08:58 PM
It splits into 3 useful functions which are much easier to understand:

Code: [Select]
'12:07 08/04/2022
'POINT INSIDE POLY

function inrange(float p,a,b) as int
====================================
if b<=p and p<a
  return -1
elseif a<=p and p<b
  return -1
endif
end function


function interpolate (float y,x1,x2,y1,y2) as float
===================================================
if y1=y2 'horizontal line
  y2=y1*+1.00001 'heuristic to avoid infinity issues
endif
return x1 + (x2 - x1) * (y - y1) / (y2 - y1)
end function


function pnpoly(int npol, float *xp,*yp, x, y) as int
=====================================================
indexbase 0
int i, j,c=0
for i=1 to i<npol
  j=i-1
  if inrange(y,yp[i],yp[j])
    if x < interpolate( y, xp[i], xp[j], yp[i], yp[j] )
      c= not c 'toggle
    endif
  endif
next
return c
end function

'TEST
uses console
float px={0,1,1,0}
float py={0,0,1,1}
float x,y
int n=4
def printi print "%1:   " cr : %1
def printd print "%1:   " %1 cr
printi ( float x=0.5, y=0.5)
printd ( pnpoly(n,px,py,x,y)     )
printd ( pnpoly(n,px,py,x+0.4,y) )
printd ( pnpoly(n,px,py,x,y+0.4) )
printd ( pnpoly(n,px,py,x+0.6,y) )
printd ( pnpoly(n,px,py,x,y+0.6) )
pause
Title: Re: Point in polygon
Post by: Nicola_Piano on April 27, 2022, 03:45:49 PM
Hi Charles,
I tried the function, it seems to be fine, but doing a practical application with gps coordinates unfortunately is not good ... in my opinion the original algorithm was already wrong ...
See the attachment.
Hello
Title: Re: Point in polygon
Post by: Charles Pegge on April 29, 2022, 07:49:19 PM
I recognize the algorithm as it is used to determine bounded areas for hatching or shading in CAD, but I think it is incomplete.

https://www.tutorialspoint.com/Check-if-a-given-point-lies-inside-a-Polygon
Title: Re: Point in polygon
Post by: Chris Chancellor on April 30, 2022, 08:55:16 AM
Good this is interesting!
hope you guys can sort out this algo as it is fairly useful to have a routine like this ?
Title: Re: Point in polygon
Post by: Charles Pegge on May 01, 2022, 11:08:50 AM
I have a solution using intersections in the inc/glo2/geoplanar.inc library. It is long but robust and tolerant of marginal cases. I'm still perfecting it.

The image below is a 4 sided shape with an inner shape. It is bombarded with 1000 points. Those which land 'inside' are highlighted yellow.



Title: Re: Point in polygon
Post by: Charles Pegge on May 01, 2022, 12:33:53 PM
Using very similar techniques, we can do cross-hatching on the interior.

the shape is 'scanned' by each hatching line for intersection points. These are collected for each line scan, sorted into ascending order of x, then pairs of points are used to draw the line segments in the interior.

Title: Re: Point in polygon
Post by: Nicola_Piano on May 02, 2022, 04:05:54 PM
Hi Charles,
it seems to me really fantastic what you managed to do. I saw the GeoPlanar.inc file, what is missing is an explanation of the various functions and input variables. Could you post an example?
Thanks.
Title: Re: Point in polygon
Post by: Charles Pegge on May 02, 2022, 07:00:10 PM
Hi Nicola,

I need another day or two, since I am still working on geoplanar, and maybe a few more examples. I also wanted to demonstrate Delaunay triangles, and their complement Voronoi diagrams, but that may be too ambitious in a short space of time.
Title: Re: Point in polygon
Post by: Nicola_Piano on May 03, 2022, 03:28:34 PM
Charles,
take your time ... the subject is quite peculiar and certainly needs special attention. :)
Title: Re: Point in polygon
Post by: Nicola_Piano on May 04, 2022, 05:56:59 PM
Hi Charles,
I managed with the intersected function to find the points that are internal or external to the polygon.
I counted the total intersections with all sides and I verified that they were odd, if they were even the point is external ....
it seems to work ....
I try it a little.
The intersected algorithm is exceptional. Thanks.
Title: Re: Point in polygon
Post by: Charles Pegge on May 06, 2022, 09:24:25 AM
Hi Nicola,

I'm still chasing a few anomalies like the missing hatch line:
Title: Re: Point in polygon
Post by: Nicola_Piano on May 10, 2022, 04:42:52 PM
Hi Charles,
how is your research going? I am eager to see progress.
Title: Re: Point in polygon
Post by: Charles Pegge on May 10, 2022, 11:01:31 PM
Yes, there is progress, though this 2d vector stuff is surprisingly tricky. I am mor accustomed to 3d graphics :)
Title: Re: Point in polygon
Post by: Charles Pegge on May 11, 2022, 01:02:26 PM
This is the code for InsideOoutsidePoints, your original quest. It should still work with the original Geoplanar.inc.

Code: [Select]
  % Title "Intersections Demo"
 '% Animated
 '% ScaleUp
 '% PlaceCentral
 '% AnchorCentral
  uses consoleG
  uses GLO2\GeoPlanar

  sub boundaries()
  ================
  indexbase 0

  string tab="    "
  point  p
  int    a,i,j,k

  seed=0x123567

  line  d[8]
  float f[16]
  f={-1,-1, 1,-1,  1,1, -1,1}
  for i=0 to 7
    f[i]+=rnd()*0.4 'turbulate
  next
  '
  d[0]={ f[0],f[1],f[2],f[3]}
  d[1]={ f[2],f[3],f[4],f[5]}
  d[2]={ f[4],f[5],f[6],f[7]}
  d[3]={ f[6],f[7],f[0],f[1]}
  '
  for i=0 to 7
    f[i]=f[i]*0.5 'smaller inner form
  next
  '
  d[4]={ f[0],f[1],f[2],f[3]}
  d[5]={ f[2],f[3],f[4],f[5]}
  d[6]={ f[4],f[5],f[6],f[7]}
  d[7]={ f[6],f[7],f[0],f[1]}



  'shading
  flat 'default
  '
  color 0,1,1
  'printl "Inside " tab p.x tab p.y

  scale 10,10

  move 1.5,-1
  '
  color 0,1,1
  thickness 3
  '
  for i=0 to 7
    drawline @d[i]
  next
  '
  int inside
  point pio
  line dd
  '
  float x,y
  for j=1 to 1000
    x=2*rnd
    y=2*rnd
    pio={x,y}
    gosub checkInside
  next
  exit sub
  '
  CheckInside:
  ============
  color 1,1,0
  pointsize 4
  dd={-1000,pio.y,1000,pio.y}
  inside=0
  for i=0 to 7
   a=intersected d[i],dd,p
   if a
     'drawline @dd
     'drawpoint @p
     if p.x<pio.x
       inside=1-inside 'toggle
     endif
   endif
  next
  '
  if inside
    color 1,1,0
  else
    color 0.8,0,0
  endif
  pointsize 4
  drawpoint @pio
  ret 
 



  end sub


  sub main()
  ==========
  string s
  cls 0,0,0
  pushstate
  scale 2,2
  'typeface=1
  printl "Inside / Outside Points"
  popstate
  '
  move 0,-2
  '
  PushState
  move 1,-1
  boundaries()
  PopState
  '
 
  end sub

  EndScript
Title: Re: Point in polygon
Post by: Nicola_Piano on May 12, 2022, 04:55:51 PM
Hi Charles,
did not come out the same result as yours ... :-(
Title: Re: Point in polygon
Post by: Charles Pegge on May 12, 2022, 05:23:37 PM
.Ok, Sorry about that. Here is my current Geoplanar.inc which belongs in /inc/GLO2

I think it will be quite open-ended for some time.

You can resize ConsoleG apps and take a jpeg-snapshot with Ctrl-P.
Title: Re: Point in polygon
Post by: Nicola_Piano on May 13, 2022, 03:52:21 PM
Hi Charles,
OK now.
I am attaching what I had done which also seems to work quite well, even with the old geoplanar.inc

Code: [Select]
' % Title "Console Demo"
 '% Animated
 '% ScaleUp
 '% PlaceCentral
 '% AnchorCentral
  uses consoleG
  uses console
  uses GLO2\GeoPlanar

  sub drawintersections()
  =======================
  string tab="    "
  line   d1,d2
  point  p
  sys    a
  int n,j,i,k,np,xt
  n=10 'number of poly
  np=9 'number of point to evaluate

 'points of poly
float px={41.84076349611504,41.837170604812826,41.83233185912741,41.82782814339409,41.82006567695157,41.81968232089131,41.81450678946148,41.809282912371984,41.81431509508004,41.820880300959615}
float py={12.466150927636479,12.480682808848268,12.48042560741089,12.478946699145974,12.481390112801055,12.468530040932215,12.461135499607634,12.45496266511059,12.44788962558273,12.438437472759132}

 'points to prove
float fx={41.834844,41.8300796,41.829849,41.832556,41.8284089,41.828539,41.8320845,41.8393343,41.820684}
float fy={12.471219,12.4668936,12.467641,12.463344,12.4608684,12.471123,12.4649883,12.4637847,12.479901}

for j=1 to np
xt=0 'conta il numero di intersenzioni
for i=1 to n
k=i+1
if i=n then k=1
d1={px[i],py[i],px[k],py[k]}
'd2={41.834844,12.471219,41.834844,0}
d2={fx[j],fy[j],fx[j],0}
a=intersected d1,d2,p
if a then
printl "CROSS" tab j "," i ")" tab a tab p.x tab p.y
xt++
end if
p.x=0
p.y=0
next i
if frac(xt/2)=0 then printl fx[j] "," fy[j] tab "> OUT"
if frac(xt/2)<>0 then printl fx[j] "," fy[j] tab "> IN"
next j
printl

end sub


sub main
=============================
  string s
  printl "Intersection Points"
  DrawIntersections
  pause
  exit
end sub

  EndScript

Title: Re: Point in polygon
Post by: Theo Gottwald on May 15, 2022, 07:26:13 PM
This stuff reminds me about neural networks where you can train the network to make classification of elements :-).
The advantage is that you do not need to find the formula, because the network will find the formula.
Neural network optimizations would be an interesting topic for oxygen possibly.
Because the basic underlying formulas are simple and can easily be optimzed in ASM.


Neural Networks for Classification (https://www.youtube.com/results?search_query=neural+network+classification)
Title: Re: Point in polygon
Post by: Charles Pegge on May 16, 2022, 03:53:09 AM
Hi Theo,

Delaunay triangles connect points so that no triangle may intersect another triangle. It can be solved by creating all possible lines between the points, then sorting into ascending order of length and testing each for intersections, giving priority to the shorter lines.

This uses random points, but it could be more interesting with other arrangements .
Title: Re: Point in polygon
Post by: Nicola_Piano on May 16, 2022, 12:57:46 PM
Charles,
this topic is interesting. I have to see to deepen it ...
Cheers
Title: Re: Point in polygon
Post by: Johan Klassen on May 16, 2022, 05:03:16 PM
Good day everyone :)
this is a very deep and interesting subject, you find a lot of info in the web like https://www.cs.cmu.edu/~quake/triangle.html