## Point In Polygon algorithm

Show off your games, demos and other (playable) creations.
devdev
Prole
Posts: 5
Joined: Thu Jun 15, 2017 5:52 am

### Point In Polygon algorithm

Greetings forumers! I ported from JS to Love2D this algorithm:
https://en.wikipedia.org/wiki/Point_in_polygon
(original code in JS: https://github.com/arkenidar/point_in_polygon)

Code: Select all

function love.load()
-- from: http://notebook.kulchenko.com/zerobrane/love2d-debugging
if arg[#arg] == "-debug" then require("mobdebug").start() end
end

local polygon={ {x=50,y=50}, {x=50,y=100}, {x=70,y=70}, {x=100,y=100}, {x=100,y=50}, }

function draw_polygon(polygon)
for x=0,500 do
for y=0,500 do
if in_concave_polygon(x, y, polygon) then
love.graphics.points({ {x,y} })
end
end
end
end

function set_polygon(new)
polygon = new
end

function love.draw()
draw_polygon(polygon)
end

function sign(x)
if x == 0 then
return 0
elseif x > 0 then
return 1
else
return -1
end
end

function side(px, p1, p2)
return sign((p2.x - p1.x) * (px.y - p1.y) - (p2.y - p1.y) * (px.x - p1.x))
end

function ring_index(index, size)
if index <= size then
if index < 1 then
return index+size
else
return index
end
else
return index-size
end
end

-- x,y: a given point, polygon: a polygon (list of points)
function in_convex_polygon(x, y, polygon)

-- minimum of 3 vertices
-- minimo 3 vertici oppure considerato punto esterno.
if #polygon < 3 then return false end

for i = 1, #polygon do
local i1,i2
i1 = ring_index(i, #polygon)
i2 = ring_index(i + 1, #polygon)
if side({x=x,y=y}, polygon[i1], polygon[i2]) > 0 then
return false
end
end

return true

end

function in_concave_polygon(x, y, polygon)

-- minimum of 3 vertices
-- minimo 3 vertici oppure considerato punto esterno.
if #polygon < 3 then return false end

while true do

-- find the first concavity vertex if there is one
-- trova la prima concavità se c'è, ovvero
-- il primo vertice che rende questo poligono un poligono concavo.
local p = polygon
local first_concavity_vertex = nil

for i=1,#polygon do
local i1,i2,i3
i1 = ring_index(i + 1, #polygon)
i2 = ring_index(i - 1, #polygon)
i3 = ring_index(i, #polygon)
local is_concavity = side(p[i1], p[i2], p[i3]) == 1
if is_concavity then
first_concavity_vertex = i
break
end
end

-- if no concavity is found then the polygon is convex
-- senza concavità il poligono è convesso.
-- (la gestione dei poligoni convessi è semplice
-- e si combina con questa casistica dei poligoni anche concavi
-- se almeno una concavità è presente).
if first_concavity_vertex == nil then
return in_convex_polygon(x, y, polygon)
end

-- if polygon is concave:
-- se il poligono è concavo:

-- form a "concavity triangle"
-- 1 - forma un triangolo di concavità
-- con il vertice concavo ed i 2 vertici adiacenti
local vertex_index = first_concavity_vertex

local i1, i2, i3
i1 = ring_index(vertex_index + 1, #polygon)
i2 = ring_index(vertex_index, #polygon)
i3 = ring_index(vertex_index - 1, #polygon)
local concavity = { polygon[i1], polygon[i2], polygon[i3] }

-- if point in "concavity triangle" return false (outside)
-- 2 - se il punto è nel triangolo di concavità (passo 1)
-- il punto è fuori
if in_convex_polygon(x, y, concavity) then
return false
end

-- form a new (!) polygon with the concavity removed
-- 3 - forma un nuovo poligono con questa concavità rimossa
-- (o meglio con rimosso il punto della prima concavità)

local new_polygon = {}
for i=1,#polygon do
if i ~= vertex_index then
table.insert( new_polygon, polygon[i] )
end
end

-- go to the cycle beginning using the polygon
-- with the concavity removed as a polygon to test
-- 4 - ripeti ma usando il nuovo poligono con
-- quel punto di concavità rimosso (passo 3) come poligono da testare
polygon = new_polygon

end

end

function love.keypressed(key, u)
--Debug
if key == "rctrl" then
--set to whatever key you want to use
--currently bound to RightControl key
debug.debug()
end
end

Thanks. Enjoy.
pgimeno
Party member
Posts: 3274
Joined: Sun Oct 18, 2015 2:58 pm

### Re: Point In Polygon algorithm

ivan
Party member
Posts: 1817
Joined: Fri Mar 07, 2008 1:39 pm
Contact:

### Re: Point In Polygon algorithm

side({x=x, y=y} means a new table instance for every step of your loop which is not good.

Check out my version stolen from vrld: https://2dengine.com/?p=polygons#Point_in_polygon
devdev
Prole
Posts: 5
Joined: Thu Jun 15, 2017 5:52 am

### Re: Point In Polygon algorithm

I am happy for the feedback! Interesting pointer links.
For the performance work can be done for sure. So I searched for "lua timing and performance profiling" since I am kind of new to this.
The code observation is right, so I patched this way:

Code: Select all

  local point = {x=x,y=y} -- only 1 time and not in the for-loop
for i = 1, #polygon do
local i1,i2
i1 = ring_index(i, #polygon)
i2 = ring_index(i + 1, #polygon)
if side(point, polygon[i1], polygon[i2]) > 0 then
return false
end
end

I can say this community is nice and welcoming!
darkfrei
Party member
Posts: 797
Joined: Sat Feb 08, 2020 11:09 pm

### Re: Point In Polygon algorithm

devdev wrote: Wed Mar 17, 2021 6:03 pm I am happy for the feedback! Interesting pointer links.
For the performance work can be done for sure. So I searched for "lua timing and performance profiling" since I am kind of new to this.
The code observation is right, so I patched this way:

Code: Select all

  local point = {x=x,y=y} -- only 1 time and not in the for-loop
for i = 1, #polygon do
local i1,i2
i1 = ring_index(i, #polygon)
i2 = ring_index(i + 1, #polygon)
if side(point, polygon[i1], polygon[i2]) > 0 then
return false
end
end

I can say this community is nice and welcoming!
Maybe step = 2?

Code: Select all

for i = 1, #polygon, 2 do
ivan
Party member
Posts: 1817
Joined: Fri Mar 07, 2008 1:39 pm
Contact:

### Re: Point In Polygon algorithm

If you want to use tables for each point then write:

Code: Select all

function in_convex_polygon(point, polygon)
You don't want to be converting these for no reason.

I don't like the ring_index function either.
A simpler and faster way is:

Code: Select all

local last = polygon[#polygon]
for i = 1, #polygon do
local current = polygon[i]
if side(point, last, current) > 0 then
return false
end
last = current
end
Lua is slightly different from JavaScript but the principles are the same.
Keep reading and studying other people's code and you'll get the hang of it.
monolifed
Party member
Posts: 188
Joined: Sat Feb 06, 2016 9:42 pm

### Re: Point In Polygon algorithm

pgimeno's code after further simplification

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) ~= (y2 > y) and (x - x2) < ((y - y2) * (x1 - x2)) / (y1 - y2) then
wn = wn + 1
end
end
return wn % 2 ~= 0 -- even/odd rule
end

quick explanation:
we use a horizontal ray (x + t, y) where t > 0

Code: Select all

(y1 > y) ~= (y2 > y)
checks whether y is between y1, y2

Code: Select all

(x - x2) < ((y - y2) * (x1 - x2)) / (y1 - y2)
checks whether slopes agree for some point (x + a, y), a > 0

Code: Select all

(x + a - x2) / (y - y2) = (x1 - x2) / (y1 - y2) so
(x - x2) < (x + a - x2) = ((y - y2) *  (x1 - x2)) / (y1 - y2)

Last edited by monolifed on Thu Mar 18, 2021 3:03 pm, edited 3 times in total.
pgimeno
Party member
Posts: 3274
Joined: Sun Oct 18, 2015 2:58 pm

### Re: Point In Polygon algorithm

Maybe that should have been posted in the other thread, as this one is about devdev's code, so apologies to devdev.

I crafted my code to have a few properties that I felt were important; specifically, to avoid divisions (no need of checking for NaNs, hopefully more speed), and to calculate the winding number, in order to allow other criteria for belonging to the poly than the odd/even rule. Specifically, winding number ~= 0 is also common; by this rule, overlaps are never considered outside.

Unfortunately, these two properties are lost in your simplification.
monolifed
Party member
Posts: 188
Joined: Sat Feb 06, 2016 9:42 pm

### Re: Point In Polygon algorithm

There is no possibility of NaN, (y1 > y) ~= (y2 > y) is false for y1 = y2

Your other points are right
Last edited by monolifed on Thu Mar 18, 2021 2:44 pm, edited 2 times in total.
pgimeno
Party member
Posts: 3274
Joined: Sun Oct 18, 2015 2:58 pm

### Re: Point In Polygon algorithm

monolifed wrote: Thu Mar 18, 2021 1:39 pm There is no possibility of NaN, (y1 > y) ~= (y2 > y) is false for y1 = y2
Actually there is: if the division overflows, it will return infinity. If the result is multiplied by zero, you get NaN.

Code: Select all

local x1, x2, y1, y2, y

y = 0
y1 = 1e-307
y2 = 0

x1 = 0
x2 = 100

assert((y1 > y) ~= (y2 > y))
print((y - y2) * ((x1 - x2) / (y1 - y2))) -- nan

In your case, removing the parentheses around the fraction would get rid of the possibility of NaN. But anyway, I was referring to my code, not yours. Multiplications, unlike divisions, need huge numbers in order to overflow.

### Who is online

Users browsing this forum: No registered users and 8 guests