Point in concave polygon

Code Sharing
1

I found myself in need of a point-in-polygon function, compatible with concave polygons (polygons containing internal angles greater than 90°, see http://en.wikipedia.org/wiki/Convex_polygon). My take on this can be seen below.

I used cross products to determine if a point is on the left or right side of the edges of a polygon. To make it work with concave polygons, a second polygon is defined that is ‘subtracted’ from the first. If a point is on the substracted polygon, it is defined as outside the polygon as a whole.

I can think of a function that determines these subtracted polygons automatically. But I was wondering if anyone has a more efficient take on this, maybe without using multiple polygons at all.

Thanks in advance!

2 
--# Main
function setup()
    -- Vertices are defined in a counter-clockwise order
    local add = {   vec2(0,0),
                    vec2(300,100),
                    vec2(400,400),
                    vec2(200,400),
                    vec2(50,300)   }
    
    local sub = {   vec2(300,100),
                    vec2(400,400),
                    vec2(200,300),
                    vec2(150,150)   }
    
    poly = ObjPolygon(add,sub)
    
    -- p Tracks the last touch
    p = vec2(100,100)
    parameter.watch("InOrOut")
end


function touched(touch)
    p = vec2(touch.x, touch.y)
end


function draw()
    background(255)
    
    poly:draw()
    fill(255, 0, 0, 255)
    ellipseMode(CENTER)
    ellipse(p.x,p.y, 10)
    
    -- Checks whether p is located in or out of the polygon
    InOrOut = poly:inPoly(p)
end


function table.copy(t)
    local copy = {}
    
    for i,v in ipairs(t) do
        copy[i] = v
    end
    
    for k,v in pairs(t) do
        copy.k = v
    end
    
    return copy
end
--# Polygon
Polygon = class()

function Polygon:init(vertices)
    self.vertices = vertices
    self.m = mesh()
    self.m.vertices = triangulate(self.vertices)
end


function Polygon:draw()
    self.m:draw()
end

-- inPoly determines if P is located in the polygon through calculating the cross product of the point P with every edge. If this product is positive, P is right of the edge. If negative, it is on the left.
function Polygon:inPoly(P)
    local v = table.copy(self.vertices)
    table.insert(v, v[1])
    
    for i = 2,#v do
        local p = P - v[i-1]
        local q = v[i] - v[i-1]
        local cross = p.x*q.y - p.y*q.x
        
        if cross > 0 then return false end
    end
    
    return true
end

--# ObjPolygon
ObjPolygon = class()

function ObjPolygon:init(vAdd, vSub)
    self.add = Polygon(vAdd)
    if vSub ~= nil then
        self.sub = Polygon(vSub)
    end
end


function ObjPolygon:draw()
    fill(0)
    self.add:draw()
    if self.sub ~= nil then
        fill(127)
        self.sub:draw()
    end
end


function ObjPolygon:inPoly(P)
    -- If P is in sub, it's
    if self.sub ~= nil and self.sub:inPoly(P) == true then
        return false
    end
    
    if self.add:inPoly(P) == true then
        return true
    end
    
    return false
end
3

I’ve implemented this by working out the winding number. This works providing the polygon doesn’t have self-intersections.

The basic routine is to add up the angles subtended at the point as you iterate over the vertices. If the result is 0 then you’re outside the polygon, if it is +/- 2pi (or +/- 360) then you’re inside. If it is +/- pi (+/- 180) then you’re on the edge.

4

Here’s a version that I had. Just move the circle thru the triangles.


displayMode(FULLSCREEN)

function setup()
    tx,ty=200,600
    verts1={vec2(200,300),vec2(500,300),vec2(350,600),
            vec2(200,400),vec2(200,400),vec2(600,500)}
end

function draw()
    background(40, 40, 50)
    stroke(255)
    strokeWidth(2)
    j=#verts1
    for z=1,#verts1 do
        line(verts1[z].x,verts1[z].y,verts1[j].x,verts1[j].y)
        j=z
    end    
    ellipse(tx,ty,10)       
    check() 
    fill(255)   
    text(str,350,700)  
end

function check()    -- check if points tx,ty are within triangle
    j=#verts1
    c=-1
    for i=1,#verts1 do
        a1=false
        if verts1[i].y>ty then
            a1=true
        end
        a2=false
        if verts1[j].y>ty then
            a2=true
        end
        b1=verts1[j].x-verts1[i].x
        b2=ty-verts1[i].y
        b3=(verts1[j].y-verts1[i].y)
        if a1~=a2 and tx<(b1*b2/b3+verts1[i].x) then
            c = c *-1
        end
        j=i
    end 
    str="outside"
    if c==1 then
        str="inside"
    end
end

function touched(t)
    tx=t.x+100
    ty=t.y+100
end
5

Here’s my checking code:

function ProtoTile:checkTouch(t)
    local tv = vec2(t.x,t.y)
    local a = 0
    local pts = self.points
    local n = self.npoints
    for k=1,n do
        a = a + (pts[k]-tv):angleBetween(pts[k%n+1] - tv)
    end
    if math.abs(a) > 1 then
        return true
    end
    return false
end
6

Thank you both for your input! I already implemented your method yesterday, @Andrew_Stacey. I like the simplicity of the concept, and it’s very clear why it works.

I am not planning on using it with complex polygons anytime soon, so that’s a challenge for the future.