Fast 2D point-in-polygon test

Showcase your libraries, tools and other projects that help your fellow love users.
Post Reply
RNavega
Party member
Posts: 235
Joined: Sun Aug 16, 2020 1:28 pm

Fast 2D point-in-polygon test

Post by RNavega »

This is a Lua port of W. Randolph Franklin's method for quickly testing if a point is within a 2D polygon. The polygon must be defined as a sequence of connected points, and can be convex or concave as well as self-intersecting.

preview.gif
preview.gif (24.96 KiB) Viewed 11080 times
Inline main.lua:

Code: Select all

--[[
    Fast 2D point-in-polygon test. This uses the "crossings count" (AKA raycast, even-odd) rule.
    
    Lua/LÖVE port by Rafael Navega (2020).
    
    License:
        Public Domain.
]]

io.stdout:setvbuf('no')


local testPoly
local renderPolys

local inside = false


function createPickablePolygon(points)
    -- Takes in a table with a sequence of ints for the (x, y) of each point of the polygon.
    -- Example: {x1, y1, x2, y2, x3, y3, ...}
    -- Note: no need to repeat the first point at the end of the table, the testing function
    -- already handles that.
    local poly = { }
    local lastX = points[#points-1]
    local lastY = points[#points]
    for index = 1, #points-1, 2 do
        local px = points[index]
        local py = points[index+1]
        -- Only store non-horizontal edges.
        if py ~= lastY then
            local index = #poly
            poly[index+1] = px
            poly[index+2] = py
            poly[index+3] = (lastX - px) / (lastY - py)
        end
        lastX = px
        lastY = py
    end
    return poly
end


function isPointInPolygon(x, y, poly)
    -- Takes in the x and y of the point in question, and a 'poly' table created by
    -- createPickablePolygon(). Returns true if the point is within the polygon, otherwise false.
    -- Note: the coordinates of the test point and the polygon points are all assumed to be in
    -- the same space.
    
    -- Original algorithm by W. Randolph Franklin (WRF):
    -- https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html
    
    local lastPX = poly[#poly-2]
    local lastPY = poly[#poly-1]
    local inside = false
    for index = 1, #poly, 3 do
        local px = poly[index]
        local py = poly[index+1]
        local deltaX_div_deltaY = poly[index+2]
        -- 'deltaX_div_deltaY' is a precomputed optimization. The original line is:
        -- if ((py > y) ~= (lastPY > y)) and (x < (y - py) * (lastX - px) / (lastY - py) + px) then        
        if ((py > y) ~= (lastPY > y)) and (x < (y - py) * deltaX_div_deltaY + px) then
            inside = not inside
        end
        lastPX = px
        lastPY = py
    end
    return inside
end


function love.load()
    love.window.setMode(620, 439)
    -- List of points taken from a flat mesh made in Blender.
    local points = {
        10, 10,
        31, 128,
        124, 170,
        10, 350,
        149, 174,
        149, 398,
        150, 398,
        150, 174,
        237, 249,
        237, 399,
        307, 399,
        307, 360,
        378, 360,
        378, 399,
        520, 399,
        520, 219,
        399, 219,
        610, 50,
        399, 199,
        399, 13,
        313, 81,
        295, 10
    }
    -- The script to print out that list of points is this, to be run in the Console view (just
    -- make sure the mesh is filled with an Ngon face first):
    --[[
        import bmesh
        bm = bmesh.from_edit_mesh(C.object.data)
        v = next(v for v in bm.verts if v.select)
        finalLoop = v.link_loops[0].link_loop_prev
        l = v.link_loops[0]
        while l != finalLoop:
            print(('%i, %i' % l.vert.co[:2]).replace('-', ''))
            l = l.link_loop_next
    ]]
    testPoly = createPickablePolygon(points)
    renderPolys = love.math.triangulate(unpack(points))
end


local DRAW_TRANSLATE_Y = 20


function love.mousemoved(x, y)
    inside = isPointInPolygon(x, y - DRAW_TRANSLATE_Y, testPoly)
end


function love.draw()
    love.graphics.print('(Hover the polygon with the mouse cursor)', 10, 10)
    love.graphics.translate(0, DRAW_TRANSLATE_Y)
    if inside then
        love.graphics.setColor(0.4, 0.65, 1.0)
    else
        love.graphics.setColor(0.2, 0.2, 0.6)
    end
    for index = 1, #renderPolys do
        love.graphics.polygon('fill', renderPolys[index])
    end
    love.graphics.setColor(1.0, 1.0, 1.0)    
end


function love.keypressed(key)
    if key == 'escape' then love.event.quit() end
end
Last edited by RNavega on Thu Oct 29, 2020 3:55 am, edited 1 time in total.
User avatar
pgimeno
Party member
Posts: 3541
Joined: Sun Oct 18, 2015 2:58 pm

Re: Fast 2D point-in-polygon test

Post by pgimeno »

Nice function. It's very fast for repeatedly querying the same polygon. If the polygon could change often, or the extra preprocessed polygons consume too much memory, I suggest this version based on the winding number, which takes the polygon without any preprocessing, but is a bit slower ~(20%) when your algorithm's preprocessing time is not considered, i.e. when repeatedly querying the same polygon using a preprocessed table:

Code: Select all

-- By Pedro Gimeno, donated to the public domain
function isPointInPolygon(x, y, poly)
  local x1, y1, x2, y2
  local len = #poly
  x2, y2 = poly[len - 1], poly[len]
  local wn = 0
  for idx = 1, len, 2 do
    x1, y1 = x2, y2
    x2, y2 = poly[idx], poly[idx + 1]

    if y1 > y then
      if (y2 <= y) and (x1 - x) * (y2 - y) < (x2 - x) * (y1 - y) then
        wn = wn + 1
      end
    else
      if (y2 > y) and (x1 - x) * (y2 - y) > (x2 - x) * (y1 - y) then
        wn = wn - 1
      end
    end
  end
  return wn % 2 ~= 0 -- even/odd rule
end
Second edit: fix a boundary condition
Last edited by pgimeno on Thu Mar 18, 2021 12:19 pm, edited 2 times in total.
RNavega
Party member
Posts: 235
Joined: Sun Aug 16, 2020 1:28 pm

Re: Fast 2D point-in-polygon test

Post by RNavega »

Hi @pgimeno, thanks for the alternative! Do you have a license on that, in case someone drops by from a web search looking for something to use?
User avatar
pgimeno
Party member
Posts: 3541
Joined: Sun Oct 18, 2015 2:58 pm

Re: Fast 2D point-in-polygon test

Post by pgimeno »

Oh duh. Updated the post. Thanks for the reminder. It's PD.
monolifed
Party member
Posts: 188
Joined: Sat Feb 06, 2016 9:42 pm

Re: Fast 2D point-in-polygon test

Post by monolifed »

Since you use +1 and -1, wouldn't "return wn ~= 0" be sufficient.
User avatar
pgimeno
Party member
Posts: 3541
Joined: Sun Oct 18, 2015 2:58 pm

Re: Fast 2D point-in-polygon test

Post by pgimeno »

monolifed wrote: Thu Mar 18, 2021 11:40 am Since you use +1 and -1, wouldn't "return wn ~= 0" be sufficient.
wn is the winding number (the number of times the polygon turns around the given point). It allows you to choose what rule to use for considering a point to be inside the polygon; see https://en.wikipedia.org/wiki/Even-odd_rule. RNavega used the even-odd rule, so I made mine compatible with that. But if you prefer the non-zero winding rule, you can change it to wn ~= 0.
Post Reply

Who is online

Users browsing this forum: Bing [Bot] and 10 guests