An implementation of Bridson's algorithm for blue noise

bridson_algorithm.py

import numpy

def sample_annulus(r1, r2, k):
	theta = numpy.random.rand(k)
	theta *= 2 * numpy.pi
	R = numpy.sqrt(numpy.random.rand(k))
	R *= r2 - r1
	R += r1
	S = numpy.array([numpy.cos(theta), numpy.sin(theta)]).T
	S *= R[:,None]
	return S


# Compute a blue noise point cloud on the unit square with Bridson's algorithm
def bridson_algorithm(radius, retry_count = 32):
	radius_sqr = radius ** 2
	N = [          (-2, -1), (-2, 0), (-2, 1),
	     (-1, -2), (-1, -1), (-1, 0), (-1, 1), (-1, 2),
	     ( 0, -2), ( 0, -1), ( 0, 0), ( 0, 1), ( 0, 2),
			 ( 1, -2), ( 1, -1), ( 1, 0), ( 1, 1), ( 1, 2),
	               ( 2, -1), ( 2, 0), ( 2, 1)]

	# Build the occupancy grid
	grid_step = radius / numpy.sqrt(2)
	grid_size = int(numpy.ceil((1. / grid_step)))
	grid = numpy.full((grid_size, grid_size), -1, dtype = int)

	# Generate the initial sample
	sample_list = [numpy.random.rand(2)]
	i, j = (int(u) for u in numpy.floor(sample_list[0] / grid_step))
	grid[i, j] = 0

	def is_valid_point(P):
		i, j = (int(u) for u in numpy.floor(P / grid_step))
		for u, v in N:
			if i + u >= 0 and j + v >= 0 and i + u < grid_size and j + v < grid_size:
				k = grid[i + u, j + v]
				if k >= 0:
					if numpy.sum((sample_list[k] - P) ** 2) < radius_sqr:
						return False

		return True

	# Generate all the samples one by one until full coverage is achieved
	active_list = [0]
	while len(active_list) > 0:
		# Randomly pick an active point n
		n = numpy.random.choice(active_list)

		# Uniform sampling in the r, 2r annulus around the active point
		candidate_list = sample_annulus(radius, 2 * radius, retry_count) + sample_list[n]

		# Discard all points out of bound
		candidate_list = [P for P in candidate_list if numpy.amax(numpy.fabs(P - .5)) <= .5]

		# Discard all points too close from previous samples
		candidate_list = [P for P in candidate_list if is_valid_point(P)]

		# No candidates found, deactive the point n
		if len(candidate_list) == 0:
			active_list.remove(n)
		else:
			active_list.append(len(sample_list))
			i, j = (int(u) for u in numpy.floor(candidate_list[0] / grid_step))
			grid[i, j] = len(sample_list)	
			sample_list.append(candidate_list[0])

	# Job done 
	return numpy.array(sample_list)