# Speed up fast poisson disk sampling generator in Python

I wrote a generator class which implements Fast Poisson Disk Sampling algorithm in Python. I did some optimizations like, x ** 2 -> x * x, using unpacking instead of indexing, move unpacking outside of loops and precalculating of constants (like 2 * pi), but still not very pleased with results. Is it possible to speed up it even more?

import math
import random

class PoissonDiskGenerator(object):
def __init__(self, field, r, k=30):
self.field_x, self.field_y = field
self.cell_size = math.ceil(r / math.sqrt(2))
self.grid_size_x, self.grid_size_y = math.ceil(field / self.cell_size), math.ceil(field / self.cell_size)
self.samples_grid = [
[None for y in range(math.ceil(self.field_x / self.cell_size))]
for x in range(math.ceil(self.field_y / self.cell_size))
]
x = random.uniform(0, field), random.uniform(0, field)
self.points = [x]
self.active_indices = 
self.active_iter = 1
self.tries = k
self.pi2 = 2 * math.pi

def __iter__(self):
return self

def __next__(self):
if self.active_indices:
point = self.try_place_new_point()
while not point and self.active_indices:
point = self.try_place_new_point()
if not point:
raise StopIteration
return point
else:
raise StopIteration

def try_place_new_point(self):
ref_ind = random.choice(self.active_indices)
for i in range(self.tries):
point_x, point_y = self.pick_point(self.points[ref_ind])
grid_x, grid_y = math.floor(point_x / self.cell_size), math.floor(point_y / self.cell_size)
neighbor_list = self.neighbors(grid_x, grid_y)
point_ok = True
if neighbor_list:
for neighbor in neighbor_list:
nb_x, nb_y = neighbor
if (point_x - nb_x) * (point_x - nb_x) + (point_y - nb_y) * (point_y - nb_y) < self.radius * self.radius:
point_ok = False
if point_ok:
self.points.append((point_x, point_y))
self.active_indices.append(self.active_iter)
self.samples_grid[grid_x][grid_y] = self.active_iter
self.active_iter += 1
return point_x, point_y
self.active_indices.remove(ref_ind)
return None

def pick_point(self, ref_point):
ref_x, ref_y = ref_point
while True:
pick_x, pick_y = ref_x + rho * math.cos(theta), ref_y + rho * math.sin(theta)
if 0 < pick_x < self.field_x and 0 < pick_y < self.field_y:
return pick_x, pick_y

def grid_to_point(self, grid_x, grid_y):
try:
return self.samples_grid[grid_x][grid_y]
except IndexError:
return None

def neighbors(self, grid_x, grid_y):
neighbors_list = (
self.grid_to_point(grid_x, grid_y),
self.grid_to_point(grid_x, grid_y - 1),
self.grid_to_point(grid_x, grid_y + 1),
self.grid_to_point(grid_x - 1, grid_y),
self.grid_to_point(grid_x - 1, grid_y - 1),
self.grid_to_point(grid_x - 1, grid_y + 1),
self.grid_to_point(grid_x + 1, grid_y),
self.grid_to_point(grid_x + 1, grid_y - 1),
self.grid_to_point(grid_x + 1, grid_y + 1),

self.grid_to_point(grid_x + 2, grid_y + 1),
self.grid_to_point(grid_x + 2, grid_y),
self.grid_to_point(grid_x + 2, grid_y - 1),

self.grid_to_point(grid_x + 1, grid_y + 2),
self.grid_to_point(grid_x, grid_y + 2),
self.grid_to_point(grid_x - 1, grid_y + 2),

self.grid_to_point(grid_x - 2, grid_y + 1),
self.grid_to_point(grid_x - 2, grid_y),
self.grid_to_point(grid_x - 2, grid_y - 1),

self.grid_to_point(grid_x + 1, grid_y - 2),
self.grid_to_point(grid_x, grid_y - 2),
self.grid_to_point(grid_x - 1, grid_y - 2)
)
return (self.points[ngb] for ngb in neighbors_list if ngb is not None)


Profiling code:

import cProfile
import pstats

def full_gen_run():
size = (15000, 15000)
point_gen = PoissonDiskGenerator(size, 100)
while True:
try:
next(point_gen)
except StopIteration:
break
print(len(point_gen.points))

cProfile.run('full_gen_run()', 'profile_stats')
stats = pstats.Stats('profile_stats')
stats.strip_dirs()
stats.sort_stats('tottime')
stats.print_stats('poissondisk.py:')


Visualisation code:

import pyglet
import time
from pyglet.window import key
from pyge.poissondisk import PoissonDiskGenerator

class Game(pyglet.window.Window):
SPEED = 10

def __init__(self):
super(Game, self).__init__(1280, 720)
self.size_x = 20000
self.size_y = 20000

self.set_caption(pyglet.version)
self.fps_display = pyglet.window.FPSDisplay(self)
pyglet.clock.schedule_interval(self.update, 1.0 / 60)
self.batch = pyglet.graphics.Batch()
self.viewpos = (self.size_x / 2, self.size_y / 2)
self.zoom = self.size_x / self.height

self.key_state_handler = key.KeyStateHandler()
self.push_handlers(self.key_state_handler)

self.point_gen = PoissonDiskGenerator((self.size_x, self.size_y), 100)
self.start_time = None
self.generation_done = False

def update(self, _):
if not self.generation_done:
if self.start_time is None:
self.start_time = time.perf_counter()
print('Points...')
time_good = True
start_time = time.perf_counter()
while time_good:
time_good = time.perf_counter() - start_time < 0.01
try:
point = next(self.point_gen)
except StopIteration:
self.generation_done = True
end_time = time.perf_counter()
print('OK ({:.2f} ms)'.format((end_time - self.start_time) * 1000))
break

if self.key_state_handler[key.W]:
self.viewpos = (self.viewpos, self.viewpos + 10 * self.SPEED)
if self.key_state_handler[key.S]:
self.viewpos = (self.viewpos, self.viewpos - 10 * self.SPEED)
if self.key_state_handler[key.A]:
self.viewpos = (self.viewpos - 10 * self.SPEED, self.viewpos)
if self.key_state_handler[key.D]:
self.viewpos = (self.viewpos + 10 * self.SPEED, self.viewpos)
if self.key_state_handler[key.E]:
self.zoom -= 0.01 * self.SPEED
if self.zoom < 1.0:
self.zoom = 1.0
if self.key_state_handler[key.Q]:
self.zoom += 0.01 * self.SPEED

def on_draw(self):
self.clear()
pyglet.gl.glViewport(0, 0, self.width, self.height)
pyglet.gl.glMatrixMode(pyglet.gl.GL_PROJECTION)
pyglet.gl.glOrtho(self.viewpos - self.width / 2 * self.zoom, self.viewpos + self.width / 2 * self.zoom,
self.viewpos - self.height / 2 * self.zoom, self.viewpos + self.height / 2 * self.zoom,
-1, 1)
pyglet.gl.glMatrixMode(pyglet.gl.GL_MODELVIEW)
self.batch.draw()
self.fps_display.draw()

if __name__ == '__main__':
game = Game()
pyglet.app.run() • What is a typical r for testing? Can you post some test code that will exercise this thing? – Reinderien Apr 10 '19 at 23:23
• @Reinderien I test it by calling next() until StopIteration raises, and gather profiling stats using cProfile. This code generate ~14179 points with r=100 and field=(15000, 15000) in about 7 sec. UPD: added code to question – Hadwig Apr 10 '19 at 23:36
• OK; but that's profiling, not testing. What kind of tests can you run against the output to ensure that it's correct? – Reinderien Apr 10 '19 at 23:58
• @Reinderien sorry, I didn't write any tests. I just visualize an output using pyglet. – Hadwig Apr 11 '19 at 0:08
• @Reinderien added code and screenshot you requested – Hadwig Apr 11 '19 at 0:39

Is it possible to speed up it even more?

Yes. Use Numpy. It's not really worth thinking about any other micro-optimizations until you've attempted to vectorize this thing with a proper numerical library.

Here's a tutorial on how to start out vectorizing with Numpy:

https://www.oreilly.com/library/view/python-for-data/9781449323592/ch04.html

There are many others.

• Any tips on what I can vectorize in this algorithm? – Hadwig Apr 10 '19 at 23:44
• Yes; I'm writing up an example – Reinderien Apr 10 '19 at 23:44
• As I see after some searching, Fast Poisson Disk algorithm can't be vectorized, because samples cannot be generated independently; each sample depends on the positions of the other samples. I can use sample elimination algorithm, but it's a different problem, that requires research to understand will it be faster or not. Question remains opened. – Hadwig Apr 11 '19 at 4:48
• @Hadwig: Maybe have a look at this implementation, it uses at least some numpy. Would be interesting to see if it is faster/slower than yours. – Graipher Apr 11 '19 at 11:50
• @Graipher This code spend 14.464 sec, while mine 8.158 sec, with same generation parametrs. Problem here - numpy ndarray work slower than python lists when you simply need to get/set one value many times. They very fast in vector calculations, but here it does not help at all. – Hadwig Apr 12 '19 at 2:00