Polygon math

Introduction

This is a small collection of polygon-related algorithms. We'll try to tackle the most common and practical issues through simple code examples.

Polygon representations

There are several different ways to represent polygons in Lua. One option is to build the polygon as a "flat" list of numbers that can be easily "unpack"-ed. While this is memory efficient, it's a little clumsy to iterate such a polygon: 
local poly = { 10,10, 20,10, 20,20, 10,20 }
for x = 1, #poly, 2 do
  local y = x + 1
  local px, py = poly[x], poly[y]
  ...
end
The second representation stores each vertex as a separate table. This technique takes more memory, but allows us to reference and modify vertices externally: 
local poly =
{
  { x=10,y=10 },
  { x=20,y=10 },
  { x=20,y=20 },
  { x=10,y=20 },
}
local first = poly[1]
first.x = first.x + 1
Ideally, it's good to stick to one polygonal representation. In reality, we are going to encounter different representations within the same API (Love2D is guilty of this too). Sooner or later we have to convert between formats. 
-- Converts to vertex table format
function toVecs(p)
  local t = {}
  local j = 1
  for i = 1, #p, 2 do
    t[j] = { x = p[i], y = p[i + 1] }
    j = j + 1
  end
  return t
end

-- Converts to flat list of coords
function toCoords(p)
  local t = {}
  local j = 1
  for i = 1, #p do
    local v = p[i]
    t[j], t[j + 1] = v.x, v.y
    j = j + 2
  end
  return t
end

Polygon generation

Regular polygons

Regular polygons are both equiangular and equilateral so all of their angles and sides are equal. We are going to look at several ways of generating regular polygons, based on different parameters. 

-- Generates a regular polygon given a circumradius (r)
function regularPoly(n, r)
  local out = {}
  for i = 1, n do
    local a = (i - 1)/n*pi2
    out[i] = { x = math.cos(a)*r, y = math.sin(a)*r }
  end
  return out
end

-- Generates a regular polygon, given a side length (s)
function regularPoly2(n, s)
  local r = s/(2*math.sin(pi/n))
  return regular(n, r)
end

-- Generates a regular polygon given an apothem (a)
function regularPoly3(n, a)
  local r = a/(math.cos(pi/n))
  return regular(n, r)
end

Exhibit 1: Regular, six-sided polygon (hexagon)
Red: circumradius
Blue: apothem
Black: sides

Tests

Clockwise or counter-clockwise

The signed area algorithm can also be used to determine the winding order of the polygon's vertices. When the signed area is negative, the vertices are are oriented in a clockwise order. A positive area means that the vertices are counter-clockwise. Zero area means that the polygon is degenerate or infinitely thin. 

-- Finds twice the signed area of a polygon
function signedPolyArea(p)
  local s = 0
  local n = #p
  local a = p[n]
  for i = 1, n do
    local b = p[i]
    s = s + (b.x + a.x)*(b.y - a.y)
    a = b
  end
  return s
end

-- Finds the actual area of a polygon
function polyArea(p)
  local s = signedPolyArea(p)
  return math.abs(s/2)
end

-- Checks if the winding of a polygon is counter-clockwise
function isPolyCCW(p)
  return signedPolyArea(p) > 0
end

-- Reverses the vertex winding
function polyReverse(p)
  local n = #p
  for i = 1, math.floor(n/2) do
    local i2 = n - i + 1
    p[i], p[i2] = p[i2], p[i]
  end
end

Exhibit 2: Counter-clockwise vs clockwise vertex winding

Convex or concave

All polygons are either convex or concave. In a convex polygon, every interior angle is less than or equal to 180 degrees. This makes convex polygons much simpler to deal with computationally.
The following algorithm uses the signed triangle area of each vertex to determine if the polygon is convex or concave. Basically, the algorithm checks if the winding order is consistent for every three consecutive vertices. 

-- Checks if a polygon is convex
function isPolyConvex(p)
  local ccw = isPolyCCW(p)
  local a = p[2]
  local b = p[1]
  for i = #p, 1, -1 do
    local c = p[i]
    local s = (c.x - a.x)*(b.y - a.y) - (c.y - a.y)*(b.x - a.x)
    if ccw then
      s = -s
    end
    if s > 0 then
      return false
    end
    a = b
    b = c
  end
  return true
end

Exhibit 3: Convexity check

Simple or complex

Complex polygons may have one or more intersecting edges. The following is a brute force algorithm that (in the worst case) test every pair of edges once. Note that we are using the segment intersection test, described in the previous tutorial. 

-- Checks if two segments intersect
function segmentOnSegment(a, b, c, d)
  local dx1, dy1 = b.x - a.x, b.y - a.y
  local dx2, dy2 = d.x - c.x, d.y - c.y
  local dx3, dy3 = a.x - c.x, a.y - c.y
  local q = dx1*dy2 - dy1*dx2
  if q == 0 then
    return false
  end
  local t1 = (dx2*dy3 - dy2*dx3)/q
  if t1 < 0 or t1 > 1 then
    return false
  end
  local t2 = (dx1*dy3 - dy1*dx3)/q
  if t2 < 0 or t2 > 1 then
    return false
  end
  return true
end

-- Checks if a polygon is simple
function isPolySimple(p)
  local n = #p
  local finish = n - 1
  local a = p[n]
  for i = 1, n do
    local b = p[i]
    local c = p[i + 1]
    for j = i + 2, finish do
      local d = p[j]
      if segmentOnSegment(a, b, c, d) then
        return false
      end
      c = d
    end
    a = b
    finish = n
  end
  return true
end

Exhibit 4: Self-intersection test

Intersections

Point in polygon

To check if a given point lies inside a polygon, we cast a ray from the point in any random direction. Next, we count the number of times that the same ray intersects with the edges of the polygon. If the number of odd, the point is inside the polygon and if it is even, the point must be outside of the polygon. Astonishingly, the algorithm below seems to work fine with both simple and self-intersecting polygons. Credit goes to Inside Code for this simple but elegant solution. 

-- Checks if a point is inside the polygon
function pointInPoly(p, q)
  local inside = false
  local n = #p
  local a = p[n]
  for i = 1, n do
    local b = p[i]
    if (q.x < a.x) ~= (q.y < b.y) then
      if q.x < a.x + ((q.y - a.y)/(b.y - a.y)*(b.x - a.x)) then
        inside = not inside
      end
    end
    a = b
  end
  return inside
end

Exhibit 5: Point vs polygon intersection test

Operations

Triangulation

Triangulation is a method of decomposing any polygon into triangles. Given a simple polygon with N number of vertices, the triangulation will always produce N - 2 triangles. One simple technique for triangulation is called ear clipping. The algorithm looks at each vertex and its neighbors and determines if that vertex is either an "ear" or "mouth" (technically called "reflex" or "convex"). The protruding "ears" are removed or "clipped" and the process is repeated on the resulting polygon. The algorithm even works with some "weakly simple" polygons that have shared edges. 

-- Finds twice the signed area of a triangle
function signedTriArea(p1, p2, p3)
  return (p1.x - p3.x)*(p2.y - p3.y) - (p1.y - p3.y)*(p2.x - p3.x)
end

-- Checks if a point is inside a triangle
function pointInTri(p, p1, p2, p3)
  local ox, oy = p.x, p.y
  local px1, py1 = p1.x - ox, p1.y - oy
  local px2, py2 = p2.x - ox, p2.y - oy
  local px3, py3 = p3.x - ox, p3.y - oy
  local sab = px1*py2 - py1*px2 < 0
  if sab ~= (px2*py3 - py2*px3 < 0) then
    return false
  end
  return sab == (px3*py1 - py3*px1 < 0)
end

-- Checks if the vertex is an "ear" or "mouth"
local left, right = {}, {}
local function isEar(i0, i1, i2)
  if signedTriArea(i0, i1, i2) >= 0 then
    local j1 = right[i2]
    repeat
      local j0, j2 = left[j1], right[j1]
      if signedTriArea(j0, j1, j2) <= 0 then
        if pointInTri(j1, i0, i1, i2) then
          return false
        end
      end
      j1 = j2
    until j1 == i0
    return true
  end
  return false
end

-- Triangulates a counter-clockwise polygon
function triangulatePoly(p)
  if not isPolyCCW(p) then
    polyReverse(p)
  end
  for i = 1, #p do
    local v = p[i]
    left[v], right[v] = p[i - 1], p[i + 1]
  end
  local first, last = p[1], p[#p]
  left[first], right[last] = last, first
  local out = {}
  local nskip = 0
  local i1 = first
  while #p >= 3 and nskip <= #p do
    local i0, i2 = left[i1], right[i1]
    if #p > 3 and isEar(i0, i1, i2) then
      table.insert(out, { i0, i1, i2 })
      left[i2], right[i0] = i0, i2     
      left[i1], right[i1] = nil, nil
      nskip = 0
      i1 = i0
    else
      nskip = nskip + 1
      i1 = i2
    end
  end
  return out
end

Exhibit 7: Triangulation

Tracing perimeter

Imagine that you are walking along the perimeter of a polygon taking steps of equal length. The following function finds points along the perimeter of a polygon, given a specific distance. The perimeter is iterated along each side, starting from the first vertex in the polygon. Note that this function returns nil if the given distance is greater than the total length of the perimeter. 

-- Finds a point along the perimeter of a polygon/polyline
function tracePerimeter(p, t, loop)
  local n = 0
  local e = loop and #p or #p - 1
  for i = 1, e do
    local a = p[i]
    local b = p[i + 1] or p[1]
    local dx = b.x - a.x
    local dy = b.y - a.y
    local d = math.sqrt(dx*dx + dy*dy)
    if n + d > t then
      local r = (t - n)/d
      return a.x + r*dx, a.y + r*dy
    end
    n = n + d
  end
end

Exhibit 8: Tracing