It splits into 3 useful functions which are much easier to understand:
'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
