Gfrktyl

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

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.