Poisson Disc Sampling in LuaTeX – Implementation Guide

luatexpatterntikz-pgf

In this answer by @Jake to the question Filling specified area by random dots in TikZ I learned about the poisson disc sampling technique to generate nice looking random patterns.

Jake's answer provided a link to a page where the algorithm is described in detail, and some implementations (in Java, Python and Ruby) are provided.

Of course this immediately sparked my urge to do a lualatex implementation which will allow to generate this kind of patterns from tikz.

So I did it. So this is not a real question, since I have the answer (which I will post in short), but I needed to share it 🙂

But I don't want to miss the opportunity to ask about nice applications for this technique in tikz. In my answer I'll focus on pattern generation, but perhaps someone more imaginative than me can think of more awesome showcases (for example, tag cloud generation?)

Best Answer

As promised, here is my implementation. It is composed of two files:

poisson.lua contains the implementation in Lua of the algorithm described here. The "main" function is generate_poisson() which returns a Lua array of points, each point being a pair (x,y). In addition, the function poisson_points_list() is provided, which transforms the data generated by generate_poisson() into something that can be fed to TikZ's foreach loop.
poisson.sty is only a wrapper to call the Lua function with TeX syntax. It provides the macro \poissonpointslist which receives four arguments: the width and height of the area to fill, the minimum distance allowed between points, and the "number of points" generated in some step of the algorithm, which is usually 20 or 30. Smaller values cause faster execution, but can leave some regions of the filled area with less dots.

The intended way of using this is to evaluate the macro with a given set of parameters and store the result in another macro, e.g.:

\edef\mylist{\poissonpointslist{5}{5}{1}{20}}

And then this new macro can be used as part of a \foreach loop, like the following:

\foreach \x,\y in \mylist { ... }

Depending on the code in the loop, we can generate an interesting variety of designs (see "Examples" section).

The code

poisson.sty

\directlua{dofile("poisson.lua")}
\newcommand{\poissonpointslist}[4]{
    \directlua{poisson_points_list(#1,#2,#3,#4)}
}

poisson.lua

-- This is a lua implementation of the algorithm described in
-- http://devmag.org.za/2009/05/03/poisson-disk-sampling/
--
-- The structure of the algorithm is exactly the same than in 
-- the mentioned article. Its pseudo-code snippets were translated
-- to Lua.
--
-- One detail worths an explanation, though. The article uses a 2D matrix
-- called grid to store coordinates of points. In the article, it is
-- assumed that grid elements can be accesed via grid[point], being point
-- some structure with a pair of x and y integers, so grid[point] should
-- be equivalent to grid[x,y] or grid[x][y]. This grid is assumed to be
-- initially dimensioned and filled by nils.
--
-- In my implementation the grid is dynamic, and it is an associative array
-- indexed by string keys in the form grid["(x,y)"]. The function gridToString()
-- can be used to convert a Point to its string form, so the grid is indeed
-- accesed like this: grid[gridToString(p)] being p a Point with integer
-- coordinates (which in fact is found via imageToGrid, like in the article)

-- UTILITY FUNCTIONS (used in the article, without giving implementation)
-- =====================================================================

-- RandomQueue stores values and gives them in random order
RandomQueue = {}
function RandomQueue.new ()
    return {last=-1}
end

function RandomQueue.push(q, item) 
    local last = q.last + 1
    q.last = last
    q[last] = item
end

function RandomQueue.pop(q)
    if (RandomQueue.empty(q)) then 
        return nil
    end
    local index = math.random(0,q.last) 
    -- A random index is generated. The last element
    -- is swaped with this random item, and the new
    -- last item is popped.
    local last = q.last
    item = q[index]
    q[index] = q[last]
    q[last] = nil
    q.last = last -1
    return item
end

function RandomQueue.empty(q)
    return q.last==-1
end

function RandomQueue.print(q)
    -- For debugging. Not used
    local t = {}
    for i=0, q.last do
        table.insert(t, string.format("(%f,%f)", q[i].x, q[i].y))
    end
    print (string.format("RandomQueue %d elem: %s", q.last+1, table.concat(t, " ")))
end

-- Point stores a coordinate pair
Point = {}

function Point.new(x,y)
    return {x=x, y=y}
end

-- Determines if a point is inside the rendered rectangle
function inRectangle(point, width, height)
  return (point.x>0 and point.y>0 and point.x<width and point.y<height)
end

-- Converts a point to a string representation, to be used as index in the grid
function gridToString(gridPoint)
    return string.format("(%d,%d)", gridPoint.x, gridPoint.y)
end

-- Computes the distance between two points
function distance(p1, p2)
    return math.sqrt(math.pow(p2.x-p1.x,2) + math.pow(p2.y-p1.y,2))
end

-- Prints the grid. For debugging. Not used
function printGrid(grid)
    print "==========="
    for k,v in pairs(grid) do
        print (string.format("%s:  %f, %f", k, v.x, v.y))
    end
end

-- THE FUNCTIONS GIVEN IN THE ARTICLE
-- This is the lua implementation of the pseudocode in the article


function generate_poisson(width, height, min_dist, new_points_count)
    local cellSize = min_dist/math.sqrt(2)
    local grid = {} -- Point.new(math.ceil(width/cellSize), math.ceil(height/cellSize))}
    local processList = RandomQueue.new()
    local samplePoints = {};  -- Empty list

    -- Generate the first point
    local firstPoint = Point.new(math.random()*width, math.random()*height)
    -- print (string.format("newPoint: [%f, %f]", firstPoint.x, firstPoint.y))
    RandomQueue.push(processList, firstPoint)
    table.insert(samplePoints, firstPoint)
    grid[gridToString(imageToGrid(firstPoint, cellSize))] = firstPoint

    -- Generate other points from points in queue
    while (not RandomQueue.empty(processList)) do
        -- RandomQueue.print(processList)
        -- printGrid(grid)
        local point = RandomQueue.pop(processList)
        for i=0,new_points_count do
            local newPoint = generateRandomPointAround(point, min_dist)
            -- print (string.format("newPoint: [%f, %f]", newPoint.x, newPoint.y))
            -- Check the point is in the region and not too close
            -- to other points
            if inRectangle(newPoint, width, height) and 
                not inNeighbourhood(grid, newPoint, min_dist, cellSize) then
                -- In this case, the point is accepted
                RandomQueue.push(processList, newPoint)
                table.insert(samplePoints, newPoint)
                grid[gridToString(imageToGrid(newPoint, cellSize))] = newPoint;
            end
        end
    end
    return samplePoints
end

function imageToGrid(point, cellSize)
    local gridX = math.floor(point.x/cellSize)
    local gridY = math.floor(point.y/cellSize)
    return Point.new(gridX, gridY)
end

function generateRandomPointAround(point, mindist)
    local r1 = math.random()
    local r2 = math.random()
    local radius  = mindist * (r1+1)
    local angle   = 2 * math.pi * r2
    newX = point.x + radius * math.cos(angle)
    newY = point.y + radius * math.sin(angle)
    return Point.new(newX, newY)
end

function inNeighbourhood(grid, point, mindist, cellSize)
    local gridPoint = imageToGrid(point, cellSize)
    cellsAroundPoint = squareAroundPoint(grid, gridPoint, 5)
    for k,cell in pairs(cellsAroundPoint) do
        if not (cell==nil) then
            local d = distance(cell, point)
            if distance(cell, point) < mindist then
                return true
            end
        end
    end
    return false
end

-- This one is not given in the article. It returns the
-- values of several cells around the give gridPoint
-- We are using string indexes for the grid, but if we
-- try to access to a key which is not stored, lua gives
-- nil instead of an exception, so it works as expected
-- because we get nils for cells which have no dot inside
function squareAroundPoint(grid, gridPoint, n)
    local extreme = math.floor(n/2)
    local result = {}
    for i=-extreme,extreme do
        for j=-extreme,extreme do
            ii = i + gridPoint.x
            jj = j + gridPoint.y
            data = grid[gridToString(Point.new(ii,jj))]
            if data == nil then
                repr = "nil"
            else
                repr = string.format("(%f,%f)", data.x, data.y)
            end
            table.insert(result, data)
        end
    end
    return result
end


-- Initialize random seed
math.randomseed(os.time())

-- Function to generate the list of dots in a tikz's foreach compatible syntax
function poisson_points_list(width, height, mindist, add_points)
    local data = generate_poisson(width, height, mindist, add_points)
    local str = {}
    for k,v in pairs(data) do
        table.insert(str, string.format("%f/%f", v.x, v.y))
    end
    tex.print(table.concat(str, ", "))
end

Examples

1. Dense distribution. Simply draw a dot at each coordinate:

\documentclass{article}
\usepackage{tikz}
\usepackage{poisson}
\begin{document}
\edef\mylist{\poissonpointslist{5}{5}{0.1}{20}}  % very dense, very slow
\begin{tikzpicture}
     \clip (0,0) rectangle (5,5);
     \foreach \x/\y in \mylist {
       \fill (\x,\y) circle(1pt);
     }
     \draw[very thick] (0,0) rectangle(5,5);
\end{tikzpicture}
\end{document}

Dense dots

2. Sparse distribution. Draw a solid disc with random radius at each point:

\edef\mylist{\poissonpointslist{5}{5}{0.7}{10}}  % Sparse, fast
\begin{tikzpicture}
    \clip (0,0) rectangle (5,5);
    \foreach \x/\y in \mylist {
      \pgfmathsetmacro{\radius}{1+7*rnd}
      \fill (\x,\y) circle(\radius pt);
    }
    \draw[very thick] (0,0) rectangle(5,5);
\end{tikzpicture}

Sparse random discs

3. Sparse density, use a flower pic at each coordinate randomizing size and angle

\tikzset{
    flower/.pic = {
       \foreach \a in {0,60,...,350}{
         \filldraw[fill=white, draw=black!30, rotate=\a+30] (0.8,0) ellipse(0.6 and 0.3);
         \filldraw[fill=white, draw=black!30, rotate=\a] (0.8,0) ellipse(0.8 and 0.3);
        }
       \filldraw[draw=orange, top color=yellow, bottom color=orange] (0,0) circle(0.4);
    },
}
\edef\mylist{\poissonpointslist{5}{5}{0.7}{10}}  % Sparse, fast
\begin{tikzpicture}
    \clip (0,0) rectangle (5,5);
    \foreach \x/\y in \mylist {
      \pgfmathsetmacro{\size}{0.15+0.1*rnd}
      \pgfmathsetmacro{\angle}{60*rnd}
      \path (\x,\y) pic[scale=\size, rotate=\angle] {flower};
    }
    \draw[very thick] (0,0) rectangle(5,5);
\end{tikzpicture}

Poisson flowers

4. More dense, use a bubble pic at each coordinate randomizing size

\tikzset{
    bubble/.pic = {
       \fill[top color=blue!50!cyan!40!white, bottom color=blue!40!black] circle(0.1);
       \fill[white] (110:0.06) circle(0.02);
    }
}
\edef\mylist{\poissonpointslist{5}{5}{0.2}{15}}  % Dense slow
\begin{tikzpicture}
     \clip (0,0) rectangle (5,5);
     \fill[top color=blue!60!cyan, bottom color=blue!70!black] (0,0) rectangle (5,5);
     \foreach \x/\y in \mylist {
       \pgfmathsetmacro{\size}{0.5+0.5*rnd}
       \path (\x,\y) pic[scale=\size]{bubble};
     }
     \draw[very thick] (0,0) rectangle(5,5);
\end{tikzpicture}

Bubbles

5. Same density as bubbles, use line segments with random angle and color

\edef\mylist{\poissonpointslist{5}{5}{0.2}{15}}  % Dense slow
\begin{tikzpicture}
     \clip (0,0) rectangle (5,5);
     \fill[orange!10] (0,0) rectangle (5,5);
     \foreach \x/\y in \mylist {
       \pgfmathsetmacro{\shade}{50*rnd}
       \pgfmathsetmacro{\w}{0.4+0.8*rnd}
       \draw[orange!50!black!\shade!white, line width=\w pt] (\x,\y) -- +(180*rnd:3pt);
     }
     \draw[very thick] (0,0) rectangle(5,5);
\end{tikzpicture}

Nice texture

Update: Bonus

To celebrate the bounty, I'm extending the answer with some considerations about the order of the list returned by \poissonpointslist and some ideas to change it. I also made all non-globally required identifiers local, as suggested by Aditya in a comment (this forced also a change in the order in which the functions are declared in the lua file). I'm pasting the new code at the end.

Order of the list

The list of coordinates returned by \poissonpoinstslist is in the order in which the algorithm generates the points. It start at a random point in the area, and then tries to find some neighbors which do not "collide" with other already set points. The growing pattern of the mesh of points is random, but always somewhat around a nucleus, as if it were the cristallization of a dissolution.

The pattern can be made visible if we change the color at which each dot is drawn, while processing the list, as in the following example (which also prints the index inside each node):

\documentclass{article}

% This makes always the same "random" list, for comparison purposes
\directlua{ math.randomseed(1); }  
\edef\mylist{\poissonpointslist{5}{5}{0.4}{20}}

\tikzset{
    disc/.style={
        circle, draw, fill=#1, font=\tiny,
        inner sep=0pt, minimum size=5mm,
    },
    disc/.default=white,
}

\begin{tikzpicture}
    \fill[black!20] (0,0) rectangle (5,5);
    \foreach [count=\i] \x/\y in \mylist {
        \pgfmathsetmacro{\tint}{100-\i/2}
        \node[disc=orange!50!white!\tint!black] at (\x,\y) {\i};
     }
\end{tikzpicture}

The result is the following:

Result

The first coordinate generated by the algorithm is the labelled as "1" which has the brighter fill. The second one is labelled as "2" and it is a bit darker, and so on. The growing pattern can be seen as a gradient of darkness.

In some scenarios it could be preferible that the order of the list were not dependent on the order in which the points were generated by the algorith, but instead having some other ordering. For example, if we want to draw the points from left-to-right and top-to-bottom, it would be preferable that the point labelled as "85" were indeed the first one, the "26" the second one, and so on.

In the code which is pasted at the end, I provide a new function \poissonpointslistordered which tries to do this (although the ordering is not perfect because there are not clear "rows" due to the random nature of the coordinates). Using this order, and the same code than above (only changing \poissonpointslist to \poissonpointslistordered), the result is now:

Result ordered

This ordered list can be useful to create some effects like the following:

\begin{tikzpicture}
    \fill[black!20] (0,0) rectangle (5,5);
    \foreach [count=\i] \x/\y in \mylist {
        \pgfmathsetmacro{\radius}{3+int(\i/5)/4}
        \fill[black!60] (\x,\y) circle (\radius pt);
     }
\end{tikzpicture}

Precipitated

The new code

poisson.sty

\directlua{dofile("poisson.lua")}
\newcommand{\poissonpointslist}[4]{
    \directlua{poisson_points_list(#1,#2,#3,#4)}
}
\newcommand{\poissonpointslistordered}[4]{
    \directlua{poisson_points_list_ordered(#1,#2,#3,#4)}
}

poisson.lua

-- This is a lua implementation of the algorithm described in
-- http://devmag.org.za/2009/05/03/poisson-disk-sampling/
--
-- The structure of the algorithm is exactly the same than in 
-- the mentioned article. Its pseudo-code snippets were translated
-- to Lua.
--
-- One detail worths an explanation, though. The article uses a 2D matrix
-- called grid to store coordinates of points. In the article, it is
-- assumed that grid elements can be accesed via grid[point], being point
-- some structure with a pair of x and y integers, so grid[point] should
-- be equivalent to grid[x,y] or grid[x][y]. This grid is assumed to be
-- initially dimensioned and filled by nils.
--
-- In my implementation the grid is dynamic, and it is an associative array
-- indexed by string keys in the form grid["(x,y)"]. The function gridToString()
-- can be used to convert a Point to its string form, so the grid is indeed
-- accesed like this: grid[gridToString(p)] being p a Point with integer
-- coordinates (which in fact is found via imageToGrid, like in the article)

-- UTILITY FUNCTIONS (used in the article, without giving implementation)
-- =====================================================================

-- RandomQueue stores values and gives them in random order
local RandomQueue = {}
function RandomQueue.new ()
    return {last=-1}
end

function RandomQueue.push(q, item) 
    local last = q.last + 1
    q.last = last
    q[last] = item
end

function RandomQueue.pop(q)
    if (RandomQueue.empty(q)) then 
        return nil
    end
    local index = math.random(0,q.last) 
    -- A random index is generated. The last element
    -- is swaped with this random item, and the new
    -- last item is popped.
    local last = q.last
    item = q[index]
    q[index] = q[last]
    q[last] = nil
    q.last = last -1
    return item
end

function RandomQueue.empty(q)
    return q.last==-1
end

function RandomQueue.print(q)
    -- For debugging. Not used
    local t = {}
    for i=0, q.last do
        table.insert(t, string.format("(%f,%f)", q[i].x, q[i].y))
    end
    print (string.format("RandomQueue %d elem: %s", q.last+1, table.concat(t, " ")))
end

-- Point stores a coordinate pair
local Point = {}

function Point.new(x,y)
    return {x=x, y=y}
end

-- Determines if a point is inside the rendered rectangle
local function inRectangle(point, width, height)
  return (point.x>0 and point.y>0 and point.x<width and point.y<height)
end

-- Converts a point to a string representation, to be used as index in the grid
local function gridToString(gridPoint)
    return string.format("(%d,%d)", gridPoint.x, gridPoint.y)
end

-- Computes the distance between two points
local function distance(p1, p2)
    return math.sqrt(math.pow(p2.x-p1.x,2) + math.pow(p2.y-p1.y,2))
end

-- Prints the grid. For debugging. Not used
local function printGrid(grid)
    print "==========="
    for k,v in pairs(grid) do
        print (string.format("%s:  %f, %f", k, v.x, v.y))
    end
end

-- THE FUNCTIONS GIVEN IN THE ARTICLE
local function imageToGrid(point, cellSize)
    local gridX = math.floor(point.x/cellSize)
    local gridY = math.floor(point.y/cellSize)
    return Point.new(gridX, gridY)
end

local function generateRandomPointAround(point, mindist)
    local r1 = math.random()
    local r2 = math.random()
    local radius  = mindist * (r1+1)
    local angle   = 2 * math.pi * r2
    newX = point.x + radius * math.cos(angle)
    newY = point.y + radius * math.sin(angle)
    return Point.new(newX, newY)
end

-- This one is not given in the article. It returns the
-- values of several cells around the give gridPoint
-- We are using string indexes for the grid, but if we
-- try to access to a key which is not stored, lua gives
-- nil instead of an exception, so it works as expected
-- because we get nils for cells which have no dot inside
local function squareAroundPoint(grid, gridPoint, n)
    local extreme = math.floor(n/2)
    local result = {}
    for i=-extreme,extreme do
        for j=-extreme,extreme do
            ii = i + gridPoint.x
            jj = j + gridPoint.y
            data = grid[gridToString(Point.new(ii,jj))]
            if data == nil then
                repr = "nil"
            else
                repr = string.format("(%f,%f)", data.x, data.y)
            end
            table.insert(result, data)
        end
    end
    return result
end

local function inNeighbourhood(grid, point, mindist, cellSize)
    local gridPoint = imageToGrid(point, cellSize)
    cellsAroundPoint = squareAroundPoint(grid, gridPoint, 5)
    for k,cell in pairs(cellsAroundPoint) do
        if not (cell==nil) then
            local d = distance(cell, point)
            if distance(cell, point) < mindist then
                return true
            end
        end
    end
    return false
end

-- This is the lua implementation of the pseudocode in the article
function generate_poisson(width, height, min_dist, new_points_count)
    local cellSize = min_dist/math.sqrt(2)
    local grid = {} -- Point.new(math.ceil(width/cellSize), math.ceil(height/cellSize))}
    local processList = RandomQueue.new()
    local samplePoints = {};  -- Empty list

    -- Generate the first point
    local firstPoint = Point.new(math.random()*width, math.random()*height)
    -- print (string.format("newPoint: [%f, %f]", firstPoint.x, firstPoint.y))
    RandomQueue.push(processList, firstPoint)
    table.insert(samplePoints, firstPoint)
    grid[gridToString(imageToGrid(firstPoint, cellSize))] = firstPoint

    -- Generate other points from points in queue
    while (not RandomQueue.empty(processList)) do
        -- RandomQueue.print(processList)
        -- printGrid(grid)
        local point = RandomQueue.pop(processList)
        for i=0,new_points_count do
            local newPoint = generateRandomPointAround(point, min_dist)
            -- print (string.format("newPoint: [%f, %f]", newPoint.x, newPoint.y))
            -- Check the point is in the region and not too close
            -- to other points
            if inRectangle(newPoint, width, height) and 
                not inNeighbourhood(grid, newPoint, min_dist, cellSize) then
                -- In this case, the point is accepted
                RandomQueue.push(processList, newPoint)
                table.insert(samplePoints, newPoint)
                grid[gridToString(imageToGrid(newPoint, cellSize))] = newPoint;
            end
        end
    end
    return samplePoints
end


-- Initialize random seed
math.randomseed(os.time())


-- Function to generate the list of dots in a tikz's foreach compatible syntax
function poisson_points_list(width, height, mindist, add_points)
    local data = generate_poisson(width, height, mindist, add_points)
    local str = {}
    for k,v in ipairs(data) do
        table.insert(str, string.format("%f/%f", v.x, v.y))
    end
    tex.print(table.concat(str, ", "))
end

-- Function similar to the above, but the returned list is "ordered"
-- so that the points are more or less in the left-to-right, 
-- top-to-down order
function poisson_points_list_ordered(width, height, mindist, add_points)
    local cellSize = mindist/math.sqrt(2)
    local function compare_coords(a,b)
        aa = imageToGrid(a, cellSize);
        bb = imageToGrid(b, cellSize);
        if (aa.y == bb.y) then  -- If they are around the same row
            return a.x<b.x;     -- the x coord orders them
        else                    -- if not
            return a.y>b.y;     -- the y coord orders them
        end
    end

    local data = generate_poisson(width, height, mindist, add_points)
    table.sort(data, compare_coords);
    local str = {}
    for k,v in ipairs(data) do
        table.insert(str, string.format("%f/%f", v.x, v.y))
    end
    tex.print(table.concat(str, ", "))
end