Poisson disc sampling implemented in luatex

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?)

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}{
\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
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} 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 {
}
\draw[very thick] (0,0) rectangle(5,5);
\end{tikzpicture} 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} 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} 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{\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} 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: 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: 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 {
}
\end{tikzpicture} The new code

poisson.sty

\directlua{dofile("poisson.lua")}
\newcommand{\poissonpointslist}{
\directlua{poisson_points_list(#1,#2,#3,#4)}
}
\newcommand{\poissonpointslistordered}{
\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
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
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

• Very nice! A small suggestion. You can make all the lua functions local except poisson_points_list and perhaps generate_poisson. This will prevent polluting the global name space (an important consideration once Luatex packages become popular). – Aditya Jun 19 '14 at 23:50
• BTW, you code also works perfectly with ConTeXt (after a minor change in the .sty file). Perfect. – Aditya Jun 19 '14 at 23:54
• FYI, there is a nice visualization (and discussion) of the algorithm here – Scott H. Jul 17 '14 at 17:42
• @Jake Thanks for the bounty! To celebrate it I updated the answer with a new feature. – JLDiaz Jul 18 '14 at 10:05
• @ScottH. That's one of the most wonderful webpages I've seen. – Manuel Jul 18 '14 at 10:23