Fast 2D point-in-polygon test

[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 (24.96 KiB) Viewed 15628 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

[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

[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?

[pgimeno]

Oh duh. Updated the post. Thanks for the reminder. It's PD.

[monolifed]

Since you use +1 and -1, wouldn't "return wn ~= 0" be sufficient.

[pgimeno]

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.

[maneco2]

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

I have this code bellow for check if have npc is inside polygon...work nice, but possible improve this code? performance?

Code: Select all

function InsideZone(x,y,poly)
	if (x ~= nil) and (y ~= nil) and (poly ~= nil) then
		n = #poly;
		inside = false;
		p1x,p1y = poly[1].X, poly[1].Y
		p2x,p2y = (poly[n].X),(poly[n].Y)
		if y > math.min(p1y,p2y) then
			if y <= math.max(p1y,p2y) then
				if x <= math.max(p1x,p2x) then
					if p1y ~= p2y then
						xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
					end;
					if p1x == p2x or x <= xinters then
						inside = not inside
					end;
				end;
			end;
		end;
		for i=1,n-1 do
			p2x,p2y = (poly[i+1].X),(poly[i+1].Y)
			if y > math.min(p1y,p2y) then
				if y <= math.max(p1y,p2y) then
					if x <= math.max(p1x,p2x) then
						if p1y ~= p2y then
							xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
						end;
						if p1x == p2x or x <= xinters then
							inside = not inside
						end;
					end;
				end;
			end;
			p1x,p1y = p2x,p2y
		end;
		return inside;
	end;
end;

[dusoft]

Yes, this:
https://love2d.org/wiki/Shape:testPoint

[RNavega]

Yes, this:
https://love2d.org/wiki/Shape:testPoint

That's cool dusoft, nice tip. That will rely on the fast C++ Box2D implementation.
The testing code seems to be this:
https://github.com/love2d/love/blob/mai ... e.cpp#L250
They're using a facingness test based on the polygon normals: if the point in question is on the negative side of all polygon normals, then the point is inside the polygon.

The only downside to this method is that it only works on convex polygons, not concave. Also, the Löve wiki mentions that a polygon shape can have at most 8 points. Probably a Box2D restriction.
Considering that this physics method should be very fast, if your polygons fit in these criteria (convex + 8 points at most), then that solution is the best choice. Edit: actually the Lua method in the OP is still the fastest, check out this comment: viewtopic.php?p=259108#p259108

I was also thinking, for some tests where the shape is too complex for polygons, you could have it as a bitmap mask, like a big byte array. If you use a uint8_t array from LuaJIT, you can pack 8 points in a single byte (edit: at least 8 points, I don't think any platforms that Löve supports have bytes with less than 8 bits), leading to a very cost effective mask that is very fast to sample.
You'd use something like this on a point-and-click adventure game, so you'd able to paint in the clickable zones like in AGS.

[dusoft]

You can have a body with multiple chained polygons (fixtures) so you could test iterating over all of them. But yes, concave only, the whole Box2D model mostly depends on that assumption.