Fast 2D point-in-polygon test

[RNavega]

Ehhh, here's something weird.
From my tests, that Lua-based method in the OP is slightly faster than the Box2D one. Point for LuaJIT I guess, it's probably optimizing that function amazingly. I honestly thought the C++ one would be faster.

point-in-polygon_benchmark.png (5.06 KiB) Viewed 893 times

The code I'm using is this (to disable the testing, set DO_BENCHMARK in love.load() to false):

Code: Select all

io.stdout:setvbuf('no')


local testPoly
local renderPolys
local testShape

local inside = false


function createPickablePolygon(points)
    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]

        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)
    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]
        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()
    local points = {}
    local TOTAL_POINTS = 8
    local angleStep = math.pi * 2.0 / TOTAL_POINTS
    for x = 0, TOTAL_POINTS - 1 do
        local index = x * 2 + 1
        points[index]     = math.cos(x * angleStep) * 100.0 + 200.0
        points[index + 1] = math.sin(x * angleStep) * 100.0 + 200.0
    end
    testPoly = createPickablePolygon(points)
    renderPolys = love.math.triangulate(unpack(points))
    testShape = love.physics.newPolygonShape(unpack(points))

    -- Set to false to disable the benchmark.
    local DO_BENCHMARK = true

    if DO_BENCHMARK then
        -- Sleep a bit to try to get the purest results.
        love.timer.sleep(2.0)
        local luaDelta = math.huge
        local b2Delta = math.huge
        for r = 1, 10 do
            local start = love.timer.getTime()
            for i = 1, 10000 do
                inside = isPointInPolygon(210, 210, testPoly)
            end
            local tempDelta = love.timer.getTime() - start
            if tempDelta < luaDelta then luaDelta = tempDelta end
        end
        -- Get a local reference so it's slightly faster.
        local testPoint = testShape.testPoint
        for r = 1, 10 do
            local start = love.timer.getTime()
            for i = 1, 10000 do
                inside = testPoint(testShape, 0, 0, 0, 210, 210)
            end
            local tempDelta = love.timer.getTime() - start
            if tempDelta < b2Delta then b2Delta = tempDelta end
        end
        error(string.format('Delta seconds (smaller is better):\nLua: %.8f\nBox2D: %.8f',
                            luaDelta, b2Delta))
    end
end


function love.mousemoved(x, y)
    -- Lua-based method.
    --inside = isPointInPolygon(x, y, testPoly)
    -- Box2D-based method.
    inside = testShape:testPoint(0, 0, 0, x, y)
end


function love.draw()
    love.graphics.print('Hover the polygon with the mouse cursor. Press Esc to quit.', 10, 10)
    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

Considering that the Lua method supports more flexible polygons (convex, concave, self-intersecting and with more than 8 points), if this performance result is the same in other systems than mine, then that'd make it the preferable method in general.

[pgimeno]

Ehhh, here's something weird.
From my tests, that Lua-based method in the OP is slightly faster than the Box2D one. Point for LuaJIT I guess, it's probably optimizing that function amazingly. I honestly thought the C++ one would be faster.

I'm not too surprised. In LuaJIT, external function calls via standard Lua (as opposed to via FFI) either prevent compilation, or require trace stitching (meaning, switch to interpreter mode, run the function, switch back to compiled mode and resume), both of which are slow. I already got caught by that trap with the transform functions when testing cam11.

I'm concerned about this condition in your code:

Code: Select all

        if py ~= lastY then

Does that mean that polygons with vertical horizontal segments can't be used?

It's likely that allowing the division by zero will be harmless and will make it work. Infinity will be properly propagated and the check for x < infinity should return the expected result.

I think the test can be rewritten without divisions, though. If we unfold deltaX_div_deltaY as (lastX - px) / (lastY - py), then:

x < (y - py) * (lastX - px) / (lastY - py) + px
is equivalent to
x - px < (y - py) * (lastX - px) / (lastY - py)
which is equivalent to
(x - px) * (lastY - py) < (y - py) * (lastX - px)

But beware the sign of lastY - py because the last equivalence only holds if it's positive.

Edit: In fact, that's pretty much what I did in my version.
Edit2: oops, not vertical but horizontal segments

[RNavega]

I'm not too surprised. In LuaJIT, external function calls via standard Lua (as opposed to via FFI) either prevent compilation, or require trace stitching (meaning, switch to interpreter mode, run the function, switch back to compiled mode and resume), both of which are slow. I already got caught by that trap with the transform functions when testing cam11.

I didn't know that, thanks for explaining.

I'm concerned about this condition in your code:

Code: Select all

        if py ~= lastY then

Does that mean that polygons with vertical segments can't be used?

From my tests, they can be used. For example, the polygon in the OP has some neighboring points with the same vertical coordinates (py == lastY), and the same horizontal coordinates.
In any case, after your rearranging of the terms in that comparison, it's good to know that even without that check for division-by-zero that I put there thinking it'd make it safe, the logic can still work.

I think the reason why W. Franklin's algorithm works while that check is in place is that it tests the previous and the next point, so if many colinear points are omitted, the geometry "silhouette" being tested ends up being the same.
That is, if the input points are these:

point-in-poly-01.png (2.06 KiB) Viewed 835 times

The data after that function would be something like this:

point-in-poly-02.png (1.29 KiB) Viewed 835 times

With redundant points omitted, but the deltas that are unique stored in the next points as that deltaX_div_deltaY.

[pgimeno]

I'm concerned about this condition in your code:

Code: Select all

        if py ~= lastY then

Does that mean that polygons with vertical segments can't be used?

From my tests, they can be used. For example, the polygon in the OP has some neighboring points with the same vertical coordinates (py == lastY), and the same horizontal coordinates.
In any case, after your rearranging of the terms in that comparison, it's good to know that even without that check for division-by-zero that I put there thinking it'd make it safe, the logic can still work.

I made an embarrassing mistake. I thought having the same Y meant they were vertical As you have noted, it means they're horizontal.

Yeah, since it's checking the number of horizontal crossings, the left and right edges suffice; there's no need for horizontal edges.