Voronoi for Löve

-- Voronoi library, by zwzsg

-- Usage:
-- local Voronoi=require("voronoi")-- Where "voronoi" is path and name to this file, with no .lua extension, and slashes replaced by dots
-- Voronoi:CreateAll(xsize,ysize,mindist)
-- Where xsize and ysize represents the dimension of the area to fill with voronoi
-- And mindist is the minimum distance between voronoi centers
-- If arguments are not supplied, then they use their previous value, unless it's the first time, in which case it errors out

-- Instead of CreateAll, you can do it step by step:
-- Voronoi:CreateFresh(512,512,50)
-- Voronoi:CreatePoints()
-- Voronoi:CreateLines()
-- Voronoi:CreateIntersections()
-- Voronoi:CreatePolygons()

-- If you need more than one Voronoi, I think you can do it via:
-- NewVoronoi=OldVoronoi:CreateFresh(nil,xsize,ysize,mindist)
-- But I haven't tested yet!

-- The Voronoi table, is made of two parts:
-- - The Voronoi library functions, in the hash part of the table
-- - The voronoi points, in the array part of the table
-- You can think of it as Voronoi={{x=x1,y=y1},{x=x2,y=y2},{x=x3,y=y3}...,CreateAll,...}
-- So that by using ipairs on this table you will iterate on the points
-- And yet with Voronoi:Funcname() you access the methods


-- Voronoi:CreateFresh(xsize,ysize,mindist)
-- Make it an empty table, with the methods, the sizes, the mindist, but no points


-- Voronoi:CreatePoints()
-- Seed with random center points
-- We now have something like {{x=x1,x=y1},{x=x2,y=y2},{x=x3,y=y3}...,CreateAll,...}
-- You can now access the coordinates of what will become voronoi centers with something like:
-- for _,v in ipairs(Voronoi) do
--   local x,y=v.x,v.y
--   ...
-- end


-- Voronoi:CreateLines()
-- Calculate the equation of the lines between voronois.
-- It might create too much, that is, lines between voronois that are not touching.
-- By definition, these are the lines equidistants to each pair of points.
--
-- We now have something like {{x=..,y=..},{x=..,y=..},{x=..,y=..}},
--                             Lines={{a=..,b=..,c=..,between={{...},{...}}},{a=..,b=..,c=..,between={{...},{...}}},...},
--                              ...,CreateAll,...}
-- Where .. represents numbers and ... represents any number of extra things
--
-- Each entry in Lines is a line, being represented by the a,b,c of its parametric equation
-- The "between" is a list of the two points that this line is bissecting

-- Voronoi:GetLineX(line,y)
-- This function is to help you work with the parametric representation of lines
-- It takes a line, as in, a table with {a=..,b=..,c=..}, and a number y, and outputs a number x, solution of a*x+b*y+c=0
-- Oh, and if the line is horizontal, this will divide by zero!

-- Voronoi:GetLineY(line,x)
-- This function is to help you work with the parametric representation of lines
-- It takes a line, as in, a table with {a=..,b=..,c=..}, and a number x, and outputs a number y, solution of a*x+b*y+c=0
-- Oh, and if the line is vertical, this will divide by zero!

-- Voronoi:CreateIntersections()
-- Calculate the intersection of the lines separating voronois.
-- There may be too many of them, and they aren't sorted, but some of them are actually coordinates of the Voronoi polygons!
-- Each intersection has four fields: the x, the y, then a field called Lines that is a list of the two lines it is the intersection of,
-- and a field called Centers that is a list of the three (but sometimes two) Centers it originate from.

-- Voronoi:CreatePolygons()
-- Calculate the coordinates of the polygons of each voronoi.
-- They are added to each voronoi as an array called c, that array contains a list of table, each table has a x and y.
-- So we now have something like:
-- {{['x']=..,['y']=..,['c']={{x=..,y=..},{x=..,y=..},{x=..,y=..},...},}
--  {['x']=..,['y']=..,['c']={{x=..,y=..},{x=..,y=..},{x=..,y=..},...},['Lines']={{a=..,b=..,c=..,with=..,it={...}},{a=..,b=..,c=..,with=..,it={...}},...},['Intersections']={{x=..,y=..,lines={...}},{x=..,y=..,lines={...}},...}}
--  Lines={{a=..,b=..,c=..,between={...}},{a=..,b=..,c=..,between={...}},...},
--  Intersections={{x=..,y=..,Lines={...},Centers={...}},{x=..,y=..,Lines={...},Centers={...}},...}
--  xsize=.., ysize=.., mindist=..,
--  CreateFresh=...., CreatePoints=...., CreateLines=...., CreateIntersections=...., CreatePolygons=...., CreateAll=....,  GetLineX=...., GetLineY=.... }
-- Where .. are numbers, ... are extra stuff, and .... are functions.


local function CreatePoints(Voronoi)-- Require CreateFresh to have been run
    assert(Voronoi.xsize>0 and Voronoi.ysize>0 and Voronoi.mindist>0,"Bad arguments, expected three stricly positive numbers")
    local attempts=0
    while attempts<50 do
        attempts=attempts+1
        local x=Voronoi.xsize*math.random()
        local y=Voronoi.ysize*math.random()
        local isAway=true
        for _,p in ipairs(Voronoi) do
            if (x-p.x)^2+(y-p.y)^2<Voronoi.mindist^2 then
                isAway=false
            end
        end
        if isAway then
            table.insert(Voronoi,{x=x,y=y})
            attempts=0
        end
    end
end

local function AssignTable(destination,source)
    for k,_ in pairs(destination) do
        destination[k]=nil
    end
    for k,v in pairs(source) do
        destination[k]=v
    end
    return destination
end

local function FillPolygonBoundingRect(poly)
    local xmin,ymin,xmax,ymax=math.huge,math.huge,-math.huge,-math.huge
    for i=#poly.c,1,-1 do
        xmin=math.min(xmin,poly.c[i].x)
        xmax=math.max(xmax,poly.c[i].x)
        ymin=math.min(ymin,poly.c[i].y)
        ymax=math.max(ymax,poly.c[i].y)
    end
    poly.xmin=xmin
    poly.xmax=xmax
    poly.ymin=ymin
    poly.ymax=ymax
end

-- To know if a point is inside a polygon, we count the number of intersection between
-- the segments and a ray going from the point, to the outside, let's say to the right.
-- If that number is even, it's outside, if that number is odd, it's inside.
local function IsInsidePolygon(poly,x,y)-- From http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
    if x<poly.xmin or x>poly.xmax or y<poly.ymin or x>poly.ymax then -- Speedup
        return false
    end
    local c=poly.c
    local inside=false
    local j=#c
    for i=1,#c do
        if y>=math.min(c[i].y,c[j].y) and y<math.max(c[i].y,c[j].y) and
            x<(c[j].x-c[i].x)*(y-c[i].y)/(c[j].y-c[i].y)+c[i].x then
            inside=not inside
        end
        j=i
    end
    return inside
end

local function IsInsidePointList(c,x,y)
    local inside=false
    local j=#c
    for i=1,#c do
        if y>=math.min(c[i].y,c[j].y) and y<math.max(c[i].y,c[j].y) and
            x<(c[j].x-c[i].x)*(y-c[i].y)/(c[j].y-c[i].y)+c[i].x then
            inside=not inside
        end
        j=i
    end
    return inside
end


local function GetLineY(line,x)
    -- a*x+b*y+c=0  =>  b*y=-a*x-c => y=(-a*x-c)/b
    return (-line.a*x-line.c)/line.b
end

local function GetLineX(line,y)
    -- a*x+b*y+c=0  =>  a*x=-b*y-c => x=(-b*y-c)/a
    return (-line.b*y-line.c)/line.a
end

local function GetIntersection(line1,line2)
    -- a1*x+b1*y+c1=0 and a2*x+b2*y+c2=0
    -- => a1*x+b1*y+c1 = a2*x+b2*y+c2
    -- => (a2-a1)*x+(b2-b1)*y+(c2-c1) = 0
    -- =>  (b2-b1)*y+(c2-c1) = -(a2-a1)*x
    -- =>  (b2-b1)*y+(c2-c1) = (a1-a2)*x
    -- =>  (a1-a2)*x = (b2-b1)*y+(c2-c1)
    -- =>  x = ((b2-b1)*y+(c2-c1))/(a1-a2)
    -- and then replace x by that in one the initial equation
    -- meh, too complicated, let's restart and use another way
    --
    -- a1*x+b1*y+c1=0 and a2*x+b2*y+c2=0
    -- Let's the first by b2 and the second by b1
    -- => b2*a1*x+b2*b1*y+b2*c1=0 and b1*a2*x+b1*b2*y+b1*c2=0
    -- Let's substract the second to the first
    -- => (b2*a1-b1*a2)*x+(b2*b1-b1*b2)*y+(b2*c1-b1*c2)=0
    -- => (b2*a1-b1*a2)*x=(b1*c2-b2*c1)
    -- => x=(b1*c2-b2*c1)/(b2*a1-b1*a2)
    --
    -- And for y:
    -- a1*x+b1*y+c1=0 and a2*x+b2*y+c2=0
    -- Let's the first by a2 and the second by a1
    -- => a2*a1*x+a2*b1*y+a2*c1=0 and a1*a2*x+a1*b2*y+a1*c2=0
    -- Let's substract the second to the first
    -- => (a2*a1-a1*a2)*x+(a2*b1-a1*b2)*y+(a2*c1-a1*c2)=0
    -- => (a2*b1-a1*b2)*y=(a1*c2-a2*c1)
    -- => y=(a1*c2-a2*c1)/(a2*b1-a1*b2)
    --
    -- So:
    -- x=(b1*c2-b2*c1)/(b2*a1-b1*a2)
    -- y=(a2*c1-a1*c2)/(b2*a1-b1*a2)

    local bottom=(line2.b*line1.a-line1.b*line2.a)-- The (b2*a1-b1*a2) part
    if bottom==0 then
        -- Oops! Division by zero! Let's hope it doesn't happen
        -- This is the case of parallel lines btw
        return nil
    end

    return (line1.b*line2.c-line2.b*line1.c)/bottom,(line2.a*line1.c-line1.a*line2.c)/bottom
end

local function GetInterAsTable(line1,line2)
    local bottom=(line2.b*line1.a-line1.b*line2.a)
    if bottom==0 then
        error("Parallel lines!")
    end
    return {x=(line1.b*line2.c-line2.b*line1.c)/bottom,y=(line2.a*line1.c-line1.a*line2.c)/bottom,lines={line1,line2}}
end

local function CreateLines(Voronoi)-- Requires CreateFresh to have been run!
    local xsize,ysize=Voronoi.xsize,Voronoi.ysize
    Voronoi.Lines={}
    for ka=2,#Voronoi do
        for kb=1,ka-1 do
            local pa,pb=Voronoi[ka],Voronoi[kb]
            local ax,ay,bx,by=pa.x,pa.y,pb.x,pb.y
            local xs,ys={},{}
            if bx-ax>xsize/2 then
                xs={{ax=ax,bx=bx-xsize,ha=true},{ax=ax+xsize,bx=bx,hb=true}}
            elseif bx-ax<-xsize/2 then
                xs={{ax=ax,bx=bx+xsize,ha=true},{ax=ax-xsize,bx=bx,hb=true}}
            else
                xs={{ax=ax,bx=bx,ha=true,hb=true}}
            end
            if by-ay>ysize/2 then
                ys={{ay=ay,by=by-ysize,ha=true},{ay=ay+ysize,by=by,hb=true}}
            elseif by-ay<-ysize/2 then
                ys={{ay=ay,by=by+ysize,ha=true},{ay=ay-ysize,by=by,hb=true}}
            else
                ys={{ay=ay,by=by,ha=true,hb=true}}
            end
            for _,xe in ipairs(xs) do
                local ax,bx=xe.ax,xe.bx
                for _,ye in ipairs(ys) do
                    local ay,by=ye.ay,ye.by
                    local ha=xe.ha and ye.ha
                    local hb=xe.hb and ye.hb
                    if ha or hb then
                        local a=ax-bx
                        local b=ay-by
                        local c=-a*(ax+bx)/2-b*(ay+by)/2
                        table.insert(Voronoi.Lines,{a=a,b=b,c=c,between={pa,pb},ha=ha,hb=hb})
                    end
                end
            end
        end
    end
    -- Sort the lines by order of distance between the two points they're bisecting
    table.sort(Voronoi.Lines,
        function(la,lb)
            return (la.between[2].x-la.between[1].x)^2+(la.between[2].y-la.between[1].y)^2<(lb.between[2].x-lb.between[1].x)^2+(lb.between[1].y-lb.between[1].y)^2
        end)
end

local function CreateIntersections(Voronoi)-- Requires CreateLines to have been run!
    Voronoi.Intersections={}
    for kla=2,#Voronoi.Lines do
        for klb=1,kla-1 do
            local la,lb=Voronoi.Lines[kla],Voronoi.Lines[klb]
            --local pa1,pa2=la.ha and la.between[1],la.hb and la.between[2]
            --local pb1,pb2=lb.ha and lb.between[1],lb.hb and lb.between[2]
            local pa1,pa2=la.between[1],la.between[2]
            local pb1,pb2=lb.between[1],lb.between[2]
            -- Only consider intersection if the two lines have in common one of the two points they're bisecting from
            if (pa1 and (pa1==pb1 or pa1==pb2)) or (pa2 and (pa2==pb1 or pa2==pb2)) then
                -- Discard line intersecting with its doppledanger
                if pa1~=pb1 or pa2~=pb2 then
                    local x,y=GetIntersection(la,lb)
                    if x and y then
                        local Centers={}
                        if pa1 then
                            table.insert(Centers,pa1)
                        end
                        if pa2 then
                            table.insert(Centers,pa2)
                        end
                        if pb1 and pb1~=pa1 and pb1~=pa2 then
                            table.insert(Centers,pb1)
                        end
                        if pb2 and pb2~=pa1 and pb2~=pa2 then
                            table.insert(Centers,pb2)
                        end
                        --assert(#Centers==3,"Bad Centers count "..#Centers)
                        if #Centers==3 then
                            table.insert(Voronoi.Intersections,{x=x,y=y,Lines={Voronoi.Lines[kla],Voronoi.Lines[klb]},Centers=Centers})
                        end
                    end
                end
            end
        end
    end
end


local function CreatePolygons(Voronoi)-- Requires CreateIntersections to have been run!
    local epsilon=0.001
    -- Helper function to calculate a distance in a wraparound zone
    -- the square root is left out because results are only compared between each others
    local xsize,ysize=Voronoi.xsize,Voronoi.ysize
    local function SquaredGetDistance(p1,p2)
        local dx=(p1.x-p2.x)%xsize
        local dy=(p1.y-p2.y)%ysize
        if dx>xsize/2 then
            dx=xsize-dx
        end
        if dy>ysize/2 then
            dy=ysize-dy
        end
        return dx*dx+dy*dy
    end
    -- Initialise corner list to empty for all voronois
    for _,p in ipairs(Voronoi) do
        p.c={}
    end
    -- For all intersections
    for _,it in ipairs(Voronoi.Intersections) do
        -- For all three voronoi center points that intersection is about
        assert(#it.Centers==3,"Bad center count "..#it.Centers)
        for _,pa in ipairs(it.Centers) do

            -- Both lines of the intersection must be about that voronoi center point
            local good=true
            for _,line in ipairs(it.Lines) do
                if line.between[1]~=pa and line.between[2]~=pa then
                    good=false
                end
            end
            if good then
                -- Disregard intersection happening too far from the voronoi center point
                if math.abs(it.x-pa.x)<xsize/2 and math.abs(it.y-pa.y)<ysize/2 then
                    local da=SquaredGetDistance(pa,it)
                    -- For all other voronoi center points
                    local bd=math.huge
                    for _,pb in ipairs(Voronoi) do
                        if pb~=pa then
                            local d=SquaredGetDistance(pb,it)
                            if d<bd then
                                -- Remember only the closest distance
                                bd=d
                            end
                        end
                    end
                    -- If the closest other is as close as me, then it's a corner
                    if math.abs(bd-da)<epsilon then
                        table.insert(pa.c,it)
                    end
                end
            end
        end
    end
    -- Corners have been written in disorder, order them by angle
    for _,p in ipairs(Voronoi) do
        table.sort(p.c,function(c1,c2)
                return math.atan2(c1.y-p.y,c1.x-p.x)<math.atan2(c2.y-p.y,c2.x-p.x)
            end)
    end
    -- Maybe here we should also remove double corners, in case these can happen
end

local function CreateAll(Voronoi,xsize,ysize,mindist)
    Voronoi:CreateFresh(xsize or Voronoi.xsize,ysize or Voronoi.ysize,mindist or Voronoi.mindist)
    Voronoi:CreatePoints()
    Voronoi:CreateLines()
    Voronoi:CreateIntersections()
    Voronoi:CreatePolygons()
end

-- Arguments: Self, the two dimensions of the area, and the minimum distance between points
local function CreateFresh(Voronoi,xsize,ysize,mindist)
    Voronoi=AssignTable(Voronoi or {},
        {xsize=xsize, ysize=ysize, mindist=mindist,
            CreateFresh=CreateFresh,
            CreatePoints=CreatePoints,
            CreateLines=CreateLines,
            CreateIntersections=CreateIntersections,
            CreatePolygons=CreatePolygons,
            CreateAll=CreateAll,
            GetLineX=GetLineX,
            GetLineY=GetLineY,
        })
    return Voronoi
end


return CreateFresh()