CODEWITHSUNDEEP
pythonpoissonsample-data
Adrián Briceño Aguilar
Adrián Briceño Aguilar

Reputation: 73

Generalizion 3d possion disk sampling

I want to create a Poisson disk sampling for a 3d grid. I used https://github.com/emulbreh/bridson implementation to generalize it to 3D. However it seems that I'm missing something, and I can't find the problem. This implementation follows the Robert Bridson's algorithm (https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf)

What I did was to make every 2d object into 3d and created the new point with spherical coordinates. 

from random import random
from math import cos, sin, floor, sqrt, pi, ceil


def euclidean_distance(a, b):
    dx = a[0] - b[0]
    dy = a[1] - b[1]
    dz = a[2] - b[2]

    return sqrt(dx * dx + dy * dy + dz * dz)


def poisson_disc_samples(width, height,thick, r, k=5, distance=euclidean_distance, random=random):
    tau = 2 * pi
    cellsize = r / sqrt(2)

    grid_width = int(ceil(width / cellsize))
    grid_height = int(ceil(height / cellsize))
    grid_thick = int(ceil(thick / cellsize))

    grid = [None] * (grid_width * grid_height * grid_thick)
    def grid_coords(p):
        return int(floor(p[0] / cellsize)), int(floor(p[1] / cellsize)), int(floor(p[2] / cellsize))

    def fits(p, gx, gy, gz):
        yrange = list(range(max(gy - 2, 0), min(gy + 3, grid_height)))
        zrange = list(range(max(gz - 2, 0), min(gz + 3, grid_thick)))
        for x in range(max(gx - 2, 0), min(gx + 3, grid_width)):
            for y in yrange:
                for z in zrange:
                    g = grid[x + y + z  * grid_width]
                    if g is None:
                        continue
                    if distance(p, g) <= r:
                        return False
        return True
    p = width * random(), height * random(), thick * random()
    queue = [p]
    grid_x, grid_y, grid_z = grid_coords(p)
    grid[grid_x + grid_y + grid_z * grid_width] = p

    while queue:
        qi = int(random() * len(queue))
        qx, qy, qz = queue[qi]
        queue[qi] = queue[-1]
        queue.pop()
        for _ in range(k):
            alpha = tau * random()
            theta = tau * random()
            d = r * sqrt(3 * random() + 1)
            px = qx + d * cos(alpha) * sin(theta)
            py = qy + d * sin(alpha) * sin(theta)
            pz = qz + d * cos(theta)
            if not (0 <= px < width and 0 <= py < height):
                continue
            p = (px, py, pz)
            grid_x, grid_y, grid_z = grid_coords(p)
            if not fits(p, grid_x, grid_y,grid_z):
                continue
            queue.append(p)
            grid[grid_x + grid_y + grid_z * grid_width] = p
    return [p for p in grid if p is not None]

r = 10
samples = poisson_disc_samples(100, 100, 100, r=10, random=random)

print(samples)

The results that I'm looking is just the printing of the sampled 3D points, however the program doesn't print any errors and I feel that because it doesn't find the samples and afterwards the computer crashes.

I think it might be the function fits(it is the one responsible for checking if the grids are empty or not) that is causing the problem.

Upvotes: 2

Views: 985

Answers (2)

Felix Zimmermann
Felix Zimmermann

Reputation: 392

I stumbled upon this code via a google search for a python 3d implementation of this algorighm. The posted code does indeed not work! There are at least two issues:

    grid[grid_x + grid_y + grid_z * grid_width] = p

In this line (and similar lines) it is not correctly indexing the grid, which should be a 3d array (or a list replaceing a 3d array). Easiest solution is to use an actual numpy 3d array and store indexes to a flat list of points in it (that't whats suggested in the paper)

d = r * sqrt(3 * random() + 1)
px = qx + d * cos(alpha) * sin(theta)
py = qy + d * sin(alpha) * sin(theta)
pz = qz + d * cos(theta)

this should uniform sample between r and 2r in three dimensions, but I believe those samples will be biased to smaller distances.

Upvotes: 0

Sam Mason
Sam Mason

Reputation: 16214

I think your changes to the code is fine... my only concern with the original is calling queue.pop() and not breaking out of the loop when a point has been found looks a bit suspicious. if you care about the order in which samples are emitted then this will be different to the published algorithm, but otherwise it probably doesn't matter

your main problem is just trying to generate a lot of points and the method has to do quite a bit of work per-point, so it takes a while to complete! if I change the code to yield points as they're generated (i.e. add yield p near queue.append(p)) instead of just returning them all at the end then I get what looks like reasonable output

Upvotes: 1

Related Questions