poissonDiscSampler.numpy.py

# translation to python for grasshopper
# original: https://gist.github.com/mbostock/22fd67be72552774736d

# import pdb
import numpy

def poissonDiscSampler(width, height, radius):
  k = 30 # maximum number of samples before rejection
  radius2 = radius * radius
  R = 3 * radius2
  cellSize = radius * numpy.sqrt(0.5)
  gridWidth = int(numpy.ceil(width / cellSize))
  gridHeight = int(numpy.ceil(height / cellSize))
  grid = numpy.zeros((gridWidth * gridHeight, 2))
  queue = []
  queueSize = 0
  sampleSize = 0
  def poisson():
    if not sampleSize:
      return sample(numpy.random.random() * width, numpy.random.random() * height)

    # Pick a random existing sample and remove it from the queue.
    while queueSize:
      i = numpy.floor(numpy.random.random() * queueSize)
      s = queue[i]
      # Make a new candidate between [radius, 2 * radius] from the existing sample.
      for j in range(0, k):
        a = 2 * numpy.pi * numpy.random.random()
        r = numpy.sqrt(numpy.random.random() * R + radius2)
        x = s[0] + r * numpy.cos(a)
        y = s[1] + r * numpy.sin(a)
        # Reject candidates that are outside the allowed extent,
        # or closer than 2 * radius to any existing sample.
        if 0 <= x and x < width and 0 <= y and y < height and far(x, y):
          return sample(x, y)
      queue[i] = queue[--queueSize]
      queue.length = queueSize

  def far(x, y):
    i = numpy.floor(x / cellSize)
    j = numpy.floor(y / cellSize)
    i0 = numpy.max(i - 2, 0)
    j0 = numpy.max(j - 2, 0)
    i1 = numpy.min(i + 3, gridWidth)
    j1 = numpy.min(j + 3, gridHeight)
    for j in range(j0, j1):
      o = j * gridWidth
      for i in range(i0, i1):
        s = grid[o + i];
        if s:
          dx = s[0] - x
          dy = s[1] - y
          if dx * dx + dy * dy < radius2:
            return False
    return True

  def sample(x, y):
    s = [x, y]
    queue.append(s)

    # pdb.set_trace()

    grid[gridWidth * int(numpy.floor(y / cellSize)) + int(numpy.floor(x / cellSize))] = s
    ++sampleSize
    ++queueSize
    return s

  return poisson

# execute the sampling

sample = poissonDiscSampler(width, height, cellSize)
samples = []
s = sample()

while True:
  s = sample()
  print s
  if s and len(samples) < maxPoints:
    samples.append(s)
  else:
    break

print samples

poissonDiscSampler.py

# translation to python for grasshopper
# original: https://gist.github.com/mbostock/22fd67be72552774736d

# import pdb
#import numpy
import time
import math
import random

random.seed(time.time())

def poissonDiscSampler(width, height, radius):
  k = 30 # maximum number of samples before rejection
  radius2 = radius * radius
  R = 3 * radius2
  cellSize = radius * 0.70710678118654757 # numpy.sqrt(0.5)
  gridWidth = int(math.ceil(width / cellSize))
  gridHeight = int(math.ceil(height / cellSize))
  grid = [[None, None]] * (gridWidth * gridHeight)
  queue = []
  queueSize = 0
  sampleSize = 0
  def poisson():
    if not sampleSize:
      return sample(random.random() * width, random.random() * height)

    # Pick a random existing sample and remove it from the queue.
    while queueSize:
      i = math.floor(random.random() * queueSize)
      s = queue[i]
      # Make a new candidate between [radius, 2 * radius] from the existing sample.
      for j in range(0, k):
        a = 2 * math.pi * random.random()
        r = math.sqrt(random.random() * R + radius2)
        x = s[0] + r * math.cos(a)
        y = s[1] + r * math.sin(a)
        # Reject candidates that are outside the allowed extent,
        # or closer than 2 * radius to any existing sample.
        if 0 <= x and x < width and 0 <= y and y < height and far(x, y):
          return sample(x, y)
      queue[i] = queue[--queueSize]
      queue.length = queueSize

  def far(x, y):
    i = math.floor(x / cellSize)
    j = math.floor(y / cellSize)
    i0 = max(i - 2, 0)
    j0 = max(j - 2, 0)
    i1 = min(i + 3, gridWidth)
    j1 = min(j + 3, gridHeight)
    for j in range(j0, j1):
      o = j * gridWidth
      for i in range(i0, i1):
        s = grid[o + i];
        if s:
          dx = s[0] - x
          dy = s[1] - y
          if dx * dx + dy * dy < radius2:
            return False
    return True

  def sample(x, y):
    s = [x, y]
    queue.append(s)

    # pdb.set_trace()

    grid[gridWidth * int(math.floor(y / cellSize)) + int(math.floor(x / cellSize))] = s
    ++sampleSize
    ++queueSize
    return s

  return poisson

# execute the sampling

sample = poissonDiscSampler(width, height, cellSize)
samples = []

while True:
  s = sample()
  print s
  if s and len(samples) < maxPoints:
    samples.append(s)
  else:
    break