-
Notifications
You must be signed in to change notification settings - Fork 12
/
poisson.py
145 lines (125 loc) · 6.14 KB
/
poisson.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
# https://github.com/scipython/scipython_maths/blob/master/poisson_disc_sampled_noise/poisson.py
import numpy as np
# For mathematical details of this algorithm, please see the blog
# article at https://scipython.com/blog/poisson-disc-sampling-in-python/
# Christian Hill, March 2017.
class PoissonDisc(object):
"""A class for generating two-dimensional Possion (blue) noise)."""
def __init__(self, rng, width=50, height=50, r=1, k=30):
self.rng = rng
self.width, self.height = width, height
self.r = r
self.k = k
# Cell side length
self.a = r/np.sqrt(2)
# Number of cells in the x- and y-directions of the grid
self.nx, self.ny = int(width / self.a) + 1, int(height / self.a) + 1
self.reset()
def reset(self):
"""Reset the cells dictionary."""
# A list of coordinates in the grid of cells
coords_list = [(ix, iy) for ix in range(self.nx)
for iy in range(self.ny)]
# Initilalize the dictionary of cells: each key is a cell's coordinates
# the corresponding value is the index of that cell's point's
# coordinates in the samples list (or None if the cell is empty).
self.cells = {coords: None for coords in coords_list}
def get_cell_coords(self, pt):
"""Get the coordinates of the cell that pt = (x,y) falls in."""
return int(pt[0] // self.a), int(pt[1] // self.a)
def get_neighbours(self, coords):
"""Return the indexes of points in cells neighbouring cell at coords.
For the cell at coords = (x,y), return the indexes of points in the
cells with neighbouring coordinates illustrated below: ie those cells
that could contain points closer than r.
ooo
ooooo
ooXoo
ooooo
ooo
"""
dxdy = [(-1,-2),(0,-2),(1,-2),(-2,-1),(-1,-1),(0,-1),(1,-1),(2,-1),
(-2,0),(-1,0),(1,0),(2,0),(-2,1),(-1,1),(0,1),(1,1),(2,1),
(-1,2),(0,2),(1,2),(0,0)]
neighbours = []
for dx, dy in dxdy:
neighbour_coords = coords[0] + dx, coords[1] + dy
if not (0 <= neighbour_coords[0] < self.nx and
0 <= neighbour_coords[1] < self.ny):
# We're off the grid: no neighbours here.
continue
neighbour_cell = self.cells[neighbour_coords]
if neighbour_cell is not None:
# This cell is occupied: store the index of the contained point
neighbours.append(neighbour_cell)
return neighbours
def point_valid(self, pt):
"""Is pt a valid point to emit as a sample?
It must be no closer than r from any other point: check the cells in
its immediate neighbourhood.
"""
cell_coords = self.get_cell_coords(pt)
for idx in self.get_neighbours(cell_coords):
nearby_pt = self.samples[idx]
# Squared distance between candidate point, pt, and this nearby_pt.
distance2 = (nearby_pt[0]-pt[0])**2 + (nearby_pt[1]-pt[1])**2
if distance2 < self.r**2:
# The points are too close, so pt is not a candidate.
return False
# All points tested: if we're here, pt is valid
return True
def get_point(self, refpt):
"""Try to find a candidate point near refpt to emit in the sample.
We draw up to k points from the annulus of inner radius r, outer radius
2r around the reference point, refpt. If none of them are suitable
(because they're too close to existing points in the sample), return
False. Otherwise, return the pt.
"""
i = 0
while i < self.k:
rho, theta = (self.rng.uniform(self.r, 2*self.r),
self.rng.uniform(0, 2*np.pi))
pt = refpt[0] + rho*np.cos(theta), refpt[1] + rho*np.sin(theta)
if not (0 < pt[0] < self.width and 0 < pt[1] < self.height):
# This point falls outside the domain, so try again.
continue
if self.point_valid(pt):
return pt
i += 1
# We failed to find a suitable point in the vicinity of refpt.
return False
def sample(self):
"""Poisson disc random sampling in 2D.
Draw random samples on the domain width x height such that no two
samples are closer than r apart. The parameter k determines the
maximum number of candidate points to be chosen around each reference
point before removing it from the "active" list.
"""
# Pick a random point to start with.
pt = (self.rng.uniform(0, self.width),
self.rng.uniform(0, self.height))
self.samples = [pt]
# Our first sample is indexed at 0 in the samples list...
self.cells[self.get_cell_coords(pt)] = 0
# and it is active, in the sense that we're going to look for more
# points in its neighbourhood.
active = [0]
# As long as there are points in the active list, keep looking for
# samples.
while active:
# choose a random "reference" point from the active list.
idx = self.rng.choice(active)
refpt = self.samples[idx]
# Try to pick a new point relative to the reference point.
pt = self.get_point(refpt)
if pt:
# Point pt is valid: add it to samples list and mark as active
self.samples.append(pt)
nsamples = len(self.samples) - 1
active.append(nsamples)
self.cells[self.get_cell_coords(pt)] = nsamples
else:
# We had to give up looking for valid points near refpt, so
# remove it from the list of "active" points.
active.remove(idx)
return self.samples