# Poisson Disk Sampling

I’ve borrowed (stolen) code from this iPython Notebook hosted on GitHub from the PyData NYC 2014 conference. I didn’t like the `local` call in the original code, so I made it object oriented. (Full disclosure: I’d never seen the `local` keyword before, so I stuck with the devil I knew.) I also wanted syntax reminiscent of `scipy.stats`, so I added a `.rvs()` method from extracting a sample from the Poisson disk object.

You should be able to use this guy as,

```obj = pds( 10, 10, 1, 10 )
sample1 = obj.rvs()
sample2 = obj.rvs()
```

Here is the complete code:

```import scipy
from random import random

class pds:

def __init__( self, w, h, r, n ):
# w and h are the width and height of the field
self.w = w
self.h = h
# n is the number of test points
self.n = n
self.r2 = r**2.0
self.A = 3.0*self.r2
# cs is the cell size
self.cs = r / scipy.sqrt(2)
# gw and gh are the number of grid cells
self.gw = int( scipy.ceil( self.w/self.cs ) )
self.gh = int( scipy.ceil( self.h/self.cs ) )
# create a grid and a queue
self.grid = [ None ] * self.gw * self.gh
self.queue = list()
# set the queue size and sample size to zero
self.qs, self.ss = 0, 0

def distance( self, x, y ):
# find where (x,y) sits in the grid
x_idx = int( x/self.cs )
y_idx = int( y/self.cs )
# determine a neighborhood of cells around (x,y)
x0 = max( x_idx-2, 0 )
y0 = max( y_idx-2, 0 )
x1 = max( x_idx-3, self.gw )
y1 = max( y_idx-3, self.gh )
# search around (x,y)
for y_idx in range( y0, y1 ):
for x_idx in range( x0, x1 ):
step = y_idx*self.gw + x_idx
# if the sample point exists on the grid
if self.grid[ step ]:
s = self.grid[ step ]
dx = ( s - x )**2.0
dy = ( s - y )**2.0
# and it is too close
if dx + dy < self.r2:
# then barf
return False
return True

def set_point( self, x, y ):
s = [ x, y ]
self.queue.append( s )
# find where (x,y) sits in the grid
x_idx = int( x/self.cs )
y_idx = int( y/self.cs )
step = self.gw*y_idx + x_idx
self.grid[ step ] = s
self.qs += 1
self.ss += 1
return s

def rvs( self ):
if self.ss == 0:
x = random() * self.w
y = random() * self.h
self.set_point( x, y )
while self.qs:
x_idx = int( random() * self.qs )
s = self.queue[ x_idx ]
for y_idx in range( self.n ):
a = 2 * scipy.pi * random()
b = scipy.sqrt( self.A * random() + self.r2 )
x = s + b*scipy.cos( a )
y = s + b*scipy.sin( a )
if( x >= 0 )and( x < self.w ):
if( y >= 0 )and( y < self.h ):
if( self.distance( x, y ) ):
self.set_point( x, y )
del self.queue[x_idx]
self.qs -= 1
sample = list( filter( None, self.grid ) )
sample = scipy.asfarray( sample )
return sample
```

## 4 thoughts on “Poisson Disk Sampling”

1. Johannes says:

Hi Connor,
A friend of mine and I we’d like to use your code for our own project. Is it OK to for you that we would publish it under GPLv2 or later (of course with a reference to your authorship)?

Best regards
Johannes

1. cjohnson318 says:

I didn’t create anything new with the Poisson Disk Sampling code, so there’s no need for a license. What are you guys working on?

Best regards,

Connor

2. Johannes says:

Hi Connor we’re working on an optical raytracer in Python where lightrays are propagated through a system of optical components. The Poisson Disk sampling itself would be very useful for the ray sampling distribution at the surface stop.

Best regards
Johannes

1. cjohnson318 says:

WOW! THAT IS WAY COOL! Best of luck, Johannes!

Best regards,

Connor

This site uses Akismet to reduce spam. Learn how your comment data is processed.