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

Example 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

Example 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

Example 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

Example 4: Self-intersection test

Intersections

Point in polygon

The following code snippet is from Hardon Collider (written by VRLD). Surprisingly, the code seems to work fine with both simple and self-intersecting polygons. 

-- Checks if an edge cuts the ray
function cutRay(a, b, q)
  return ((a.y > q.y and b.y < q.y) or (a.y < q.y and b.y > q.y))
    and (q.x - a.x < (q.y - a.y)*(b.x - a.x)/(b.y - a.y))
end

-- Checks if the ray crosses boundary from interior to exterior
function crossBoundary(a, b, q)
  return (a.y == q.y and a.x > q.x and b.y < q.y) or
    (b.y == q.y and b.x > q.x and a.y < q.y)
end

-- 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 cutRay(a, b, q) or crossBoundary(a, b, q) then
      inside = not inside
    end
    a = b
  end
  return inside
end

Example 5: Point vs polygon intersection test

Operations

Polygon offsetting

Polygon offsetting is an operation that contracts or expands the polygon. At first glance, offsetting may look similar to scaling, but the resulting polygon is quite different. The following algorithm offsets each edge of the source polygon, in the direction of its normal. Adjacent edges are then "joined" using the "lineVsLine" function. The code will only work with convex polygons and some concave polygon. In real life cases, the length of the mitered corners needs to be clamped especially when dealing with "spiky" polygons. 

-- Offsets the given segment (a,b) by width (w)
function segmentOffset(a, b, w)
  local dx, dy = b.x - a.x, b.y - a.y
  local d = math.sqrt(dx*dx + dy*dy)
  assert(d > 0, "degenerate segment")
  local ox, oy = dx/d*w, dy/d*w
  return { x = a.x + oy, y = a.y - ox }, { x = b.x + oy, y = b.y - ox }
end

-- Find the intersection point of two lines
function lineVsline(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 d = dx1*dy2 - dy1*dx2
  if d ~= 0 then
    local t1 = (dx2*dy3 - dy2*dx3)/d
    return { x = a[1] + t1*dx1, y = a[2] + t1*dy1 }
  end
end

-- Offsets the given polygon (p) by width (w)
local tail, head = {}, {}
function polyOffset(p, w, loop)
  assert(w ~= 0, "invalid offset width")
  -- offset segments
  local n = #p
  local q = p[1]
  for i = n, 1, -1 do
    local b = p[i]
    head[i], tail[i] = segmentOffset(q, b, w)
    q = b
  end
  -- join intersections
  local t = {}
  local a = head[n]
  local b = tail[n]
  for i = 1, n do
    local c = head[i]
    local d = tail[i]
    t[i] = lineVsline(a, b, c, d) or b
    a = c
    b = d
  end
  -- loose ends
  if not loop then
    t[1] = head[1]
    t[n] = tail[n - 1]
  end
  return t
end

Example 6:
Black: original polygon
Red: offset polygon

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 ab = px1*py2 - py1*px2
  local px3, py3 = p3.x - ox, p3.y - oy
  local bc = px2*py3 - py2*px3
  local sab = ab < 0
  if sab ~= (bc < 0) then
    return false
  end
  local ca = px3*py1 - py3*px1
  return sab == (ca < 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

Example 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

Example 8: Tracing