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

Post by devdev »

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.
User avatar
ivan
Party member
Posts: 1822
Joined: Fri Mar 07, 2008 1:39 pm
Contact:

Re: Point In Polygon algorithm

Post by ivan »

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

Post by devdev »

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! :)
User avatar
darkfrei
Party member
Posts: 814
Joined: Sat Feb 08, 2020 11:09 pm

Re: Point In Polygon algorithm

Post by darkfrei »

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
User avatar
ivan
Party member
Posts: 1822
Joined: Fri Mar 07, 2008 1:39 pm
Contact:

Re: Point In Polygon algorithm

Post by ivan »

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

Post by monolifed »

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.
User avatar
pgimeno
Party member
Posts: 3285
Joined: Sun Oct 18, 2015 2:58 pm

Re: Point In Polygon algorithm

Post by pgimeno »

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

Post by monolifed »

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.
User avatar
pgimeno
Party member
Posts: 3285
Joined: Sun Oct 18, 2015 2:58 pm

Re: Point In Polygon algorithm

Post by pgimeno »

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.
Post Reply

Who is online

Users browsing this forum: No registered users and 33 guests