2
\$\begingroup\$

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[0] / self.cell_size), math.ceil(field[1] / 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[0]), random.uniform(0, field[1])
        self.points = [x]
        self.active_indices = [0]
        self.active_iter = 1
        self.tries = k
        self.radius = r
        self.radius2 = 2 * r
        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:
            rho, theta = random.uniform(self.radius, self.radius2), random.uniform(0, self.pi2)
            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
                self.batch.add(1, pyglet.gl.GL_POINTS, None, ('v2f', point))

        if self.key_state_handler[key.W]:
            self.viewpos = (self.viewpos[0], self.viewpos[1] + 10 * self.SPEED)
        if self.key_state_handler[key.S]:
            self.viewpos = (self.viewpos[0], self.viewpos[1] - 10 * self.SPEED)
        if self.key_state_handler[key.A]:
            self.viewpos = (self.viewpos[0] - 10 * self.SPEED, self.viewpos[1])
        if self.key_state_handler[key.D]:
            self.viewpos = (self.viewpos[0] + 10 * self.SPEED, self.viewpos[1])
        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.glLoadIdentity()
        pyglet.gl.glOrtho(self.viewpos[0] - self.width / 2 * self.zoom, self.viewpos[0] + self.width / 2 * self.zoom,
                          self.viewpos[1] - self.height / 2 * self.zoom, self.viewpos[1] + 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()

screenshot of visualisation

\$\endgroup\$
  • \$\begingroup\$ What is a typical r for testing? Can you post some test code that will exercise this thing? \$\endgroup\$ – Reinderien Apr 10 at 23:23
  • \$\begingroup\$ @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 \$\endgroup\$ – Hadwig Apr 10 at 23:36
  • \$\begingroup\$ OK; but that's profiling, not testing. What kind of tests can you run against the output to ensure that it's correct? \$\endgroup\$ – Reinderien Apr 10 at 23:58
  • \$\begingroup\$ @Reinderien sorry, I didn't write any tests. I just visualize an output using pyglet. \$\endgroup\$ – Hadwig Apr 11 at 0:08
  • 1
    \$\begingroup\$ @Reinderien added code and screenshot you requested \$\endgroup\$ – Hadwig Apr 11 at 0:39
1
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Any tips on what I can vectorize in this algorithm? \$\endgroup\$ – Hadwig Apr 10 at 23:44
  • \$\begingroup\$ Yes; I'm writing up an example \$\endgroup\$ – Reinderien Apr 10 at 23:44
  • \$\begingroup\$ 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. \$\endgroup\$ – Hadwig Apr 11 at 4:48
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ – Graipher Apr 11 at 11:50
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ – Hadwig Apr 12 at 2:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.