#!/usr/bin/env lua --------------- -- ## Delaunay, Lua module for convex polygon triangulation -- @author Roland Yonaba -- @copyright 2013-2016 -- @license MIT -- @script delaunay -- ================ -- Private helpers -- ================ local setmetatable = setmetatable local tostring = tostring local assert = assert local unpack = unpack or table.unpack local remove = table.remove local sqrt = math.sqrt local max = math.max -- Internal class constructor local class = function(...) local klass = {} klass.__index = klass klass.__call = function(_,...) return klass:new(...) end function klass:new(...) local instance = setmetatable({}, klass) klass.__init(instance, ...) return instance end return setmetatable(klass,{__call = klass.__call}) end -- Triangle semi-perimeter by Heron's formula local function quatCross(a, b, c) local p = (a + b + c) * (a + b - c) * (a - b + c) * (-a + b + c) return sqrt(p) end -- Cross product (p1-p2, p2-p3) local function crossProduct(p1, p2, p3) local x1, x2 = p2.x - p1.x, p3.x - p2.x local y1, y2 = p2.y - p1.y, p3.y - p2.y return x1 * y2 - y1 * x2 end -- Checks if angle (p1-p2-p3) is flat local function isFlatAngle(p1, p2, p3) return (crossProduct(p1, p2, p3) == 0) end -- ================ -- Module classes -- ================ --- `Edge` class -- @type Edge local Edge = class() Edge.__eq = function(a, b) return (a.p1 == b.p1 and a.p2 == b.p2) end Edge.__tostring = function(e) return (('Edge :\n %s\n %s'):format(tostring(e.p1), tostring(e.p2))) end --- Creates a new `Edge` -- @name Edge:new -- @param p1 a `Point` -- @param p2 a `Point` -- @return a new `Edge` -- @usage -- local Delaunay = require 'Delaunay' -- local Edge = Delaunay.Edge -- local Point = Delaunay.Point -- local e = Edge:new(Point(1,1), Point(2,5)) -- local e = Edge(Point(1,1), Point(2,5)) -- Alias to Edge.new -- print(e) -- print the edge members p1 and p2 -- function Edge:__init(p1, p2) self.p1, self.p2 = p1, p2 end --- Test if `otherEdge` is similar to self. It does not take into account the direction. -- @param otherEdge an `Edge` -- @return `true` or `false` -- @usage -- local e1 = Edge(Point(1,1), Point(2,5)) -- local e2 = Edge(Point(2,5), Point(1,1)) -- print(e1:same(e2)) --> true -- print(e1 == e2)) --> false, == operator considers the direction -- function Edge:same(otherEdge) return ((self.p1 == otherEdge.p1) and (self.p2 == otherEdge.p2)) or ((self.p1 == otherEdge.p2) and (self.p2 == otherEdge.p1)) end --- Returns the length. -- @return the length of self -- @usage -- local e = Edge(Point(), Point(10,0)) -- print(e:length()) --> 10 -- function Edge:length() return self.p1:dist(self.p2) end --- Returns the midpoint coordinates. -- @return the x-coordinate of self midpoint -- @return the y-coordinate of self midpoint -- @usage -- local e = Edge(Point(), Point(10,0)) -- print(e:getMidPoint()) --> 5, 0 -- function Edge:getMidPoint() local x = self.p1.x + (self.p2.x - self.p1.x) / 2 local y = self.p1.y + (self.p2.y - self.p1.y) / 2 return x, y end --- Point class -- @type Point local Point = class() Point.__eq = function(a,b) return (a.x == b.x and a.y == b.y) end Point.__tostring = function(p) return ('Point (%s) x: %.2f y: %.2f'):format(p.id, p.x, p.y) end --- Creates a new `Point` -- @name Point:new -- @param x the x-coordinate -- @param y the y-coordinate -- @return a new `Point` -- @usage -- local Delaunay = require 'Delaunay' -- local Point = Delaunay.Point -- local p = Point:new(1,1) -- local p = Point(1,1) -- Alias to Point.new -- print(p) -- print the point members x and y -- function Point:__init(x, y) self.x, self.y, self.id = x or 0, y or 0, '?' end --- Returns the square distance to another `Point`. -- @param p a `Point` -- @return the square distance from self to `p`. -- @usage -- local p1, p2 = Point(), Point(1,1) -- print(p1:dist2(p2)) --> 2 -- function Point:dist2(p) local dx, dy = (self.x - p.x), (self.y - p.y) return dx * dx + dy * dy end --- Returns the distance to another `Point`. -- @param p a `Point` -- @return the distance from self to `p`. -- @usage -- local p1, p2 = Point(), Point(1,1) -- print(p1:dist2(p2)) --> 1.4142135623731 -- function Point:dist(p) return sqrt(self:dist2(p)) end --- Checks if self lies into the bounds of a circle -- @param cx the x-coordinate of the circle center -- @param cy the y-coordinate of the circle center -- @param r the radius of the circle -- @return `true` or `false` -- @usage -- local p = Point() -- print(p:isInCircle(0,0,1)) --> true -- function Point:isInCircle(cx, cy, r) local dx = (cx - self.x) local dy = (cy - self.y) return ((dx * dx + dy * dy) <= (r * r)) end --- `Triangle` class -- @type Triangle local Triangle = class() Triangle.__tostring = function(t) return (('Triangle: \n %s\n %s\n %s') :format(tostring(t.p1), tostring(t.p2), tostring(t.p3))) end --- Creates a new `Triangle` -- @name Triangle:new -- @param p1 a `Point` -- @param p2 a `Point` -- @param p3 a `Point` -- @return a new `Triangle` -- @usage -- local Delaunay = require 'Delaunay' -- local Triangle = Delaunay.Triangle -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle:new(p1, p2, p3) -- local t = Triangle(p1, p2, p3) -- Alias to Triangle.new -- print(t) -- print the triangle members p1, p2 and p3 -- function Triangle:__init(p1, p2, p3) assert(not isFlatAngle(p1, p2, p3), ("angle (p1, p2, p3) is flat:\n %s\n %s\n %s") :format(tostring(p1), tostring(p2), tostring(p3))) self.p1, self.p2, self.p3 = p1, p2, p3 self.e1, self.e2, self.e3 = Edge(p1, p2), Edge(p2, p3), Edge(p3, p1) end --- Checks if the triangle is defined clockwise (sequence p1-p2-p3) -- @return `true` or `false` -- @usage -- local p1, p2, p3 = Point(), Point(1,1), Point(2,0) -- local t = Triangle(p1, p2, p3) -- print(t:isCW()) --> true -- function Triangle:isCW() return (crossProduct(self.p1, self.p2, self.p3) < 0) end --- Checks if the triangle is defined counter-clockwise (sequence p1-p2-p3) -- @return `true` or `false` -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:isCCW()) --> true -- function Triangle:isCCW() return (crossProduct(self.p1, self.p2, self.p3) > 0) end --- Returns the length of the edges -- @return the length of the edge p1-p2 -- @return the length of the edge p2-p3 -- @return the length of the edge p3-p1 -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:getSidesLength()) --> 2 1.4142135623731 1.4142135623731 -- function Triangle:getSidesLength() return self.e1:length(), self.e2:length(), self.e3:length() end --- Returns the coordinates of the center -- @return the x-coordinate of the center -- @return the y-coordinate of the center -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:getCenter()) --> 1 0.33333333333333 -- function Triangle:getCenter() local x = (self.p1.x + self.p2.x + self.p3.x) / 3 local y = (self.p1.y + self.p2.y + self.p3.y) / 3 return x, y end --- Returns the coordinates of the circumcircle center and its radius -- @return the x-coordinate of the circumcircle center -- @return the y-coordinate of the circumcircle center -- @return the radius of the circumcircle -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:getCircumCircle()) --> 1 0 1 -- function Triangle:getCircumCircle() local x, y = self:getCircumCenter() local r = self:getCircumRadius() return x, y, r end --- Returns the coordinates of the circumcircle center -- @return the x-coordinate of the circumcircle center -- @return the y-coordinate of the circumcircle center -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:getCircumCenter()) --> 1 0 -- function Triangle:getCircumCenter() local p1, p2, p3 = self.p1, self.p2, self.p3 local D = ( p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y)) * 2 local x = (( p1.x * p1.x + p1.y * p1.y) * (p2.y - p3.y) + ( p2.x * p2.x + p2.y * p2.y) * (p3.y - p1.y) + ( p3.x * p3.x + p3.y * p3.y) * (p1.y - p2.y)) local y = (( p1.x * p1.x + p1.y * p1.y) * (p3.x - p2.x) + ( p2.x * p2.x + p2.y * p2.y) * (p1.x - p3.x) + ( p3.x * p3.x + p3.y * p3.y) * (p2.x - p1.x)) return (x / D), (y / D) end --- Returns the radius of the circumcircle -- @return the radius of the circumcircle -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:getCircumRadius()) --> 1 -- function Triangle:getCircumRadius() local a, b, c = self:getSidesLength() return ((a * b * c) / quatCross(a, b, c)) end --- Returns the area -- @return the area -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:getArea()) --> 1 -- function Triangle:getArea() local a, b, c = self:getSidesLength() return (quatCross(a, b, c) / 4) end --- Checks if a given point lies into the triangle circumcircle -- @param p a `Point` -- @return `true` or `false` -- @usage -- local p1, p2, p3 = Point(), Point(2,0), Point(1,1) -- local t = Triangle(p1, p2, p3) -- print(t:inCircumCircle(Point(1,-1))) --> true -- function Triangle:inCircumCircle(p) return p:isInCircle(self:getCircumCircle()) end --- Delaunay module -- @section public --- Delaunay module -- @table Delaunay -- @field Point reference to the `Point` class -- @field Edge reference to the `Edge` class -- @field Triangle reference to the `Triangle` class -- @field convexMultiplier multiplier heuristic for bounding triangle calculation. When small (~1) produces convex-hull, when large, produces concave hulls. Defaults to 1000. -- @field _VERSION the version of the current module local Delaunay = { Point = Point, Edge = Edge, Triangle = Triangle, convexMultiplier = 1e3, _VERSION = "0.1" } --- Triangulates a set of given vertices -- @param ... a `vargarg` list of objects of type `Point` -- @return a set of objects of type `Triangle` -- @usage -- local Delaunay = require 'Delaunay' -- local Point = Delaunay.Point -- local p1, p2, p3, p4 = Point(), Point(2,0), Point(1,1), Point(1,-1) -- local triangles = Delaunay.triangulate(p1, p2, p3, p4) -- for i = 1, #triangles do -- print(triangles[i]) -- end -- function Delaunay.triangulate(...) local vertices = {...} local nvertices = #vertices assert(nvertices > 2, "Cannot triangulate, needs more than 3 vertices") if nvertices == 3 then return {Triangle(unpack(vertices))} end local trmax = nvertices * 4 local minX, minY = vertices[1].x, vertices[1].y local maxX, maxY = minX, minY for i = 1, #vertices do local vertex = vertices[i] vertex.id = i if vertex.x < minX then minX = vertex.x end if vertex.y < minY then minY = vertex.y end if vertex.x > maxX then maxX = vertex.x end if vertex.y > maxY then maxY = vertex.y end end local convex_mult = Delaunay.convexMultiplier local dx, dy = (maxX - minX) * convex_mult, (maxY - minY) * convex_mult local deltaMax = max(dx, dy) local midx, midy = (minX + maxX) * 0.5, (minY + maxY) * 0.5 local p1 = Point(midx - 2 * deltaMax, midy - deltaMax) local p2 = Point(midx, midy + 2 * deltaMax) local p3 = Point(midx + 2 * deltaMax, midy - deltaMax) p1.id, p2.id, p3.id = nvertices + 1, nvertices + 2, nvertices + 3 vertices[p1.id] = p1 vertices[p2.id] = p2 vertices[p3.id] = p3 local triangles = {} triangles[#triangles + 1] = Triangle(vertices[nvertices + 1], vertices[nvertices + 2], vertices[nvertices + 3] ) for i = 1, nvertices do local edges = {} local ntriangles = #triangles for j = #triangles, 1, -1 do local curTriangle = triangles[j] if curTriangle:inCircumCircle(vertices[i]) then edges[#edges + 1] = curTriangle.e1 edges[#edges + 1] = curTriangle.e2 edges[#edges + 1] = curTriangle.e3 remove(triangles, j) end end for j = #edges - 1, 1, -1 do for k = #edges, j + 1, -1 do if edges[j] and edges[k] and edges[j]:same(edges[k]) then remove(edges, j) remove(edges, k-1) end end end for j = 1, #edges do local n = #triangles assert(n <= trmax, "Generated more than needed triangles") triangles[n + 1] = Triangle(edges[j].p1, edges[j].p2, vertices[i]) end end for i = #triangles, 1, -1 do local triangle = triangles[i] if (triangle.p1.id > nvertices or triangle.p2.id > nvertices or triangle.p3.id > nvertices) then remove(triangles, i) end end for _ = 1,3 do remove(vertices) end return triangles end return Delaunay