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
Point in polygon

### Author Topic: Point in polygon  (Read 800 times)

0 Members and 1 Guest are viewing this topic.

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Point in polygon
« 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

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #1 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 0int 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;  endifnextreturn c;end function'TESTfloat px={0,1,1,0}float py={0,0,1,1}float x,yint n=4float x=0.5, y=0.5print pnpoly(n,px,py,x,y)print pnpoly(n,px,py,x+0.6,y)print pnpoly(n,px,py,x,y+0.6)`
« Last Edit: April 07, 2022, 08:55:04 PM by Charles Pegge »

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #2 on: April 07, 2022, 10:13:15 PM »
Thanks Charles,
tomorrow I'll try it and tell you.

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #3 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 POLYfunction inrange(float p,a,b) as int====================================if b<=p and p<a  return -1elseif a<=p and p<b  return -1endifend functionfunction interpolate (float y,x1,x2,y1,y2) as float===================================================if y1=y2 'horizontal line  y2=y1*+1.00001 'heuristic to avoid infinity issuesendifreturn x1 + (x2 - x1) * (y - y1) / (y2 - y1)end functionfunction pnpoly(int npol, float *xp,*yp, x, y) as int=====================================================indexbase 0int i, j,c=0for 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  endifnextreturn cend function'TESTuses consolefloat px={0,1,1,0}float py={0,0,1,1}float x,yint n=4def printi print "%1:   " cr : %1def printd print "%1:   " %1 crprinti ( 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`

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #4 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

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #5 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

#### Chris Chancellor

• Sr. Member
• Posts: 345
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #6 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 ?

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #7 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.

« Last Edit: May 01, 2022, 11:12:46 AM by Charles Pegge »

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #8 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.

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #9 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.

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #10 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.

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #11 on: May 03, 2022, 03:28:34 PM »
Charles,
take your time ... the subject is quite peculiar and certainly needs special attention.

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #12 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.

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 744
• User-Rate: +30/-1
##### Re: Point in polygon
« Reply #13 on: May 06, 2022, 09:24:25 AM »
Hi Nicola,

I'm still chasing a few anomalies like the missing hatch line:

#### Nicola_Piano

• Jr. Member
• Posts: 82
• User-Rate: +0/-0
##### Re: Point in polygon
« Reply #14 on: May 10, 2022, 04:42:52 PM »
Hi Charles,
how is your research going? I am eager to see progress.