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}
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}
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{\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}
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 {
\pgfmathsetmacro{\radius}{3+int(\i/5)/4}
\fill[black!60] (\x,\y) circle (\radius pt);
}
\end{tikzpicture}
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
