Point in a polygon

[BasicScience]

Anyone have an efficient routine for the "point in a polygon" problem?

The idea is that a set of connected points defines a closed polygon. Now, for any given new test point, is it inside our outside the polygon?

This method will be used for a complex polygon with about 300 points (from tracing an object on the screen). The question is "what is the area bounded by the polygon". Seems to me the easiest solution is to first find the smallest rectangle that will bound the entire polygon, then test each pixel in the rectangle to determine if it's inside or outside the polygon.

Thanks.

[sancho3]

Ray casting method and others
And you will find a couple of the methods already coded in this thread

[dodicat]

The sub circulate is to ensure the points lie radially so avoiding crossing lines.
You can test locally with the mouse.
Any polygon need not be closed.
I.E. a square would have four vertices, no need to superimpose the last on the first.
Press spacebar to refresh, Esc to end.

Code: Select all

Type Point
    As Long x,y
End Type

Function inpolygon(p1() As Point,Byval p2 As Point) As long 'point in a polygon-winding number method
    #macro Winder(L1,L2,p)
    ((L1.x-L2.x)*(p.y-L2.y)-(p.x-L2.x)*(L1.y-L2.y))
    #endmacro
    Dim As long index,nextindex,k=Ubound(p1)+1,wn
    For n As long=1 To Ubound(p1)
        index=n Mod k:nextindex=(n+1) Mod k
        If nextindex=0 Then nextindex=1
        If p1(index).y<=p2.y Then
            If p1(nextindex).y>p2.y Andalso  Winder(p1(index),p1(nextindex),p2)>0 Then wn+=1 
        Else
            If p1(nextindex).y<=p2.y Andalso Winder(p1(index),p1(nextindex),p2)<0 Then wn-=1
        End If
    Next n
    Return wn
End Function

Sub drawpolygon(p() As Point,Byref col As Ulong,Byval im As Any Pointer=0) 
    Dim k As Long=Ubound(p)+1
    Dim As Long index,nextindex
    Dim As Long cx,cy
    For n As Long=1 To Ubound(p)
        cx+=p(n).x:cy+=p(n).y
        index=n Mod k:nextindex=(n+1) Mod k
        If nextindex=0 Then nextindex=1
        Line im,(p(index).x,p(index).y)-(p(nextindex).x,p(nextindex).y),col
    Next
    cx/=Ubound(p):cy/=Ubound(p)
    Paint (cx,cy),col,col
End Sub

Sub circulate(p() As Point)'sort the random points radially 
    #macro Circlesort() 
    '  bubblesort
    For p1 As Long  = Lbound(p) To Ubound(p)-1
        For p2 As Long  = p1 + 1 To Ubound(p)
            If Atan2(p(p1).y-c.y,p(p1).x-c.x)< Atan2(p(p2).y-c.y,p(p2).x-c.x) Then
                Swap p(p1),p(p2)
            End If
        Next p2
    Next p1
    #endmacro
    Dim As Point C '--centroid of points
    Dim As Long counter
    For n As Long=Lbound(p) To Ubound(p)
        counter+=1
        c.x+=p(n).x
        c.y+=p(n).y
    Next n
    c.x=c.x/counter
    c.y=c.y/counter
    CircleSort()
End Sub

redim as point polygon()
dim as long mx,my,area
dim as string i
screen 20,32
do
    area=0
   redim polygon(1 to 3+rnd*100)
    for n as long=1 to ubound(polygon)
        polygon(n).x=rnd*1024
        polygon(n).y=rnd*768
    next
     circulate(polygon())
    dim as ulong clr=rgb(rnd*255,rnd*255,rnd*255)
    'find pixels inside
    for x as long=0 to 1024
        for y as long=0 to 768
            if inpolygon(polygon(),type(x,y)) then area+=1
        next
        next
do
    i=inkey
    getmouse(mx,my)
    
    
    screenlock
    cls
    
   
   drawpolygon(polygon(),clr)
   locate 2,2
   print iif(inpolygon(polygon(),type(mx,my)),"INSIDE","OUTSIDE")
   print "Area in pixels = ";area
   print "Number of vertices = ";ubound(polygon)
   screenunlock
   sleep 1,1
loop until len(i)

loop until i=chr(27)

    


 

[dodicat]

Or, a very simple way using the colour of the polygon

Code: Select all

#define getpixel(_x,_y)     *cptr(ulong ptr,row + (_y)*pitch + (_x) shl 2)

Type Point
    As Long x,y
End Type



Sub drawpolygon(p() As Point,Byref col As Ulong,Byval im As Any Pointer=0) 
    Dim k As Long=Ubound(p)+1
    Dim As Long index,nextindex
    Dim As Long cx,cy
    For n As Long=1 To Ubound(p)
        cx+=p(n).x:cy+=p(n).y
        index=n Mod k:nextindex=(n+1) Mod k
        If nextindex=0 Then nextindex=1
        Line im,(p(index).x,p(index).y)-(p(nextindex).x,p(nextindex).y),col
    Next
    cx/=Ubound(p):cy/=Ubound(p)
    Paint (cx,cy),col,col
End Sub

Sub circulate(p() As Point)'sort the random points radially 
    #macro Circlesort() 
    '  bubblesort
    For p1 As Long  = Lbound(p) To Ubound(p)-1
        For p2 As Long  = p1 + 1 To Ubound(p)
            If Atan2(p(p1).y-c.y,p(p1).x-c.x)< Atan2(p(p2).y-c.y,p(p2).x-c.x) Then
                Swap p(p1),p(p2)
            End If
        Next p2
    Next p1
    #endmacro
    Dim As Point C '--centroid of points
    Dim As Long counter
    For n As Long=Lbound(p) To Ubound(p)
        counter+=1
        c.x+=p(n).x
        c.y+=p(n).y
    Next n
    c.x=c.x/counter
    c.y=c.y/counter
    CircleSort()
End Sub



screen 20,32
redim as point polygon()
dim as long mx,my,area
dim as string i

dim  as any ptr row
dim  as integer pitch
row=screenptr
Screeninfo ,,,,.pitch
do
    
   redim polygon(1 to 3+rnd*100)
    for n as long=1 to ubound(polygon)
        polygon(n).x=rnd*1024
        polygon(n).y=rnd*768
    next
     circulate(polygon())
    dim as ulong clr=rgb(rnd*255,rnd*255,rnd*255)
    
do
    area=0
    i=inkey
    getmouse(mx,my)
    
    
    screenlock
    cls
    
   
   drawpolygon(polygon(),clr)
   'find pixels inside
   for x as long=0 to 1024-1
        for y as long=0 to 768-1
            if getpixel(x,y)=clr then area+=1
        next
        next
   
   locate 2,2
   print iif(getpixel(abs(mx),abs(my))=clr,"INSIDE","OUTSIDE")
   print "Area in pixels = ";area
   print "Number of vertices = ";ubound(polygon)
   screenunlock
   sleep 1,1
loop until len(i)

loop until i=chr(27)

    
 

[paul doe]

Anyone have an efficient routine for the "point in a polygon" problem?

The idea is that a set of connected points defines a closed polygon. Now, for any given new test point, is it inside our outside the polygon?

Yes, there are 2 well known solutions: the winding number and the even-odd rule tracing. Dodicat already gave you the first one, but for completeness, here are the two of them:

Code: Select all

type point2d
	as integer x
	as integer y
end type

enum point_result
	outside
	inside
end enum

#define min( x, y )   iif( x < y, x, y )
#define max( x, y )   iif( x > y, x, y )
#define isLeft( p0, p1, p2 ) ( ( p1.x - p0.x ) * ( p2.y - p0.y ) - ( p2.x - p0.x ) * ( p1.y - p0.y ) )

'' winding number rule
function wn_pnPoly( v() as point2d, byval p as point2d ) as point_result
	dim as integer n = ubound(v)
	dim as integer wn = 0
	
	for i as integer = 0 to n-1
		if( v(i).y <= p.y ) then
			if( v( i + 1 ).y > p.y ) then
				if( isleft( v( i ), v( i + 1 ), p ) > 0 ) then
					wn += 1
				end if
			end if
		else
				if( v( i + 1 ).y <= p.y ) then
					if( isleft( v( i ), v( i + 1 ), p ) < 0 ) then
						wn -= 1
					end if
				end if
		end if
	next
	
	return( iif( wn > 0, inside, outside ) )	
end function

'' even-odd rule
function insidePoly( poly() as point2D, byval p as point2D ) as point_result
  dim as integer n
  dim as point2D p1, p2
  
  dim as integer counter
  dim as double xinters
  
  p1 = poly( 0 )
  counter = 0
  
  n = ubound( poly )
  
  for i as integer = 0 to n
    p2 = poly( i mod n )
    
    if( p.y > min( p1.y, p2.y ) ) then
      if( p.y <= max( p1.y, p2.y ) ) then
        if( p.x <= max( p1.x, p2.x ) )then
          if( p1.y <> p2.y ) then
              xinters = ( p.y - p1.y ) * ( p2.x - p1.x ) / ( p2.y - p1.y ) + p1.x
              
              if( ( p1.x = p2.x ) orElse ( p.x <= xinters ) ) then
                counter += 1
              end if
          end if
        end if
      end if
    end if
    
    p1 = p2
  next
  
  if( ( counter mod 2 ) = 0 ) then
    '' outside
    return( outside )
  else
    '' inside
    return( inside )
  end if
end function

Which one of them is more efficient depends on what you plan to do with said polygon. For a single point query, the winding number is faster. If you need to rasterize the polygon, a well-implemented even-odd rule can prove more efficient (you can perform the rasterizing and the test while you traverse the ray).

[Carlos Herrera]

The question is "what is the area bounded by the polygon". ..

So, you know coordinates of the vertices? Why not calculate the area directly?
As far as I remember, it is quite easy for not self-intersecting polygons

[BasicScience]

Many thanks for the code. All examples worked well.

At first, I did not realize there is a simple solution to the area of a polygon, given the coordinates of the verticies. This site has a great review of geometric relationships and algorithms

http://geomalgorithms.com/a01-_area.html

For those who are curious,

Area = 0.5 * SUM { x * (y[i+1] - y[i-1]) }

[dodicat]

I tried out the formula.
Looks like it is in agreement with winding number (only a fraction of a percent difference)
Thanks for sharing the formula.

Code: Select all

Type Point
    As Double x,y
End Type
'optimized a bit
Function inpolygon(p1() As Point,Byval p2 As Point) As Long
    Static As Long nxt,k
    Dim As Long wn:k=Ubound(p1)
    For n As Long=1 To k
        nxt=(n+1)
        If nxt>k Then nxt=1
        If p1(n).y<=p2.y Then
            If p1(nxt).y>p2.y Then
If ((p1(n).x-p1(nxt).x)*(p2.y-p1(nxt).y)-(p2.x-p1(nxt).x)*(p1(n).y-p1(nxt).y))>0 Then wn+=1
            End If
        Else
            If p1(nxt).y<=p2.y Then
If ((p1(n).x-p1(nxt).x)*(p2.y-p1(nxt).y)-(p2.x-p1(nxt).x)*(p1(n).y-p1(nxt).y))<0 Then wn-=1
            End If
        End If
    Next n
    Return wn
End Function


Sub drawpolygon(p() As Point,Byref col As Ulong,Byval im As Any Pointer=0)
    Dim k As Long=Ubound(p)+1
    Dim As Long index,nextindex
    Dim As Long cx,cy
    For n As Long=1 To Ubound(p)
        cx+=p(n).x:cy+=p(n).y
        index=n Mod k:nextindex=(n+1) Mod k
        If nextindex=0 Then nextindex=1
        Line im,(p(index).x,p(index).y)-(p(nextindex).x,p(nextindex).y),col
    Next
    cx/=Ubound(p):cy/=Ubound(p)
    Paint (cx,cy),col,col
End Sub

Sub circulate(p() As Point)
    Dim As Point C
    Dim As Long counter
    For n As Long=Lbound(p) To Ubound(p)
        counter+=1
        c.x+=p(n).x
        c.y+=p(n).y
    Next n
    c.x=c.x/counter
    c.y=c.y/counter
    For p1 As Long =Lbound(p) To Ubound(p)-1
        For p2 As Long =p1 + 1 To Ubound(p)
            If Atan2(p(p1).y-c.y,p(p1).x-c.x)<Atan2(p(p2).y-c.y,p(p2).x-c.x) Then
                Swap p(p1),p(p2)
            End If
        Next p2
    Next p1
End Sub

Function area2(p() As Point) As Double
    Dim As Double a
    Var j=Ubound(p)
    For i As Long=Lbound(p) To Ubound(p)
        a=a+(p(j).x+p(i).x) *(p(j).y-p(i).y)
        j=i
    Next i
    Return (a)/2
End Function

Redim As Point polygon()
Dim As Long mx,my
Dim As Double area
Dim As String i
Screen 20,32
Do
    area=0
    Redim polygon(1 To 3+Rnd*100)
    For n As Long=1 To Ubound(polygon)
        polygon(n).x=Rnd*1024
        polygon(n).y=Rnd*768
    Next
    circulate(polygon())
    Dim As Ulong clr= ((Cuint(Rnd*255) Shl 16) Or (Cuint(Rnd*255) Shl 8) Or Cuint(Rnd*255) Or &hFF000000)
    Dim As Double t=Timer
    For x As Long=0 To 1024
        For y As Long=0 To 768
            If inpolygon(polygon(),Type(x,y)) Then area+=1
        Next
    Next
    Dim As Double t2=Timer-t
    Do
        i=Inkey
        Getmouse(mx,my)
        
        
        Screenlock
        Cls
        
        
        drawpolygon(polygon(),clr)
        Locate 2,2
        Print Iif(inpolygon(polygon(),Type(mx,my)),"INSIDE","OUTSIDE")
        
        Print "Number of vertices = ";Ubound(polygon)
        Print "Area in pixels = ";area
        Print "area by formula  ";area2(polygon())
        Print "difference       ";area-area2(polygon()),((area-area2(polygon()))/area)*100;" percent"
        'Print t2
        Screenunlock
        Sleep 1,1
    Loop Until Len(i)
    
Loop Until i=Chr(27)