|
| 1 | +# SPDX-FileCopyrightText: The Threadbare Authors |
| 2 | +# SPDX-License-Identifier: MPL-2.0 |
| 3 | + |
| 4 | +class_name PoissonDiscSampler |
| 5 | +extends Object |
| 6 | + |
| 7 | +# https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf |
| 8 | +# https://www.jasondavies.com/poisson-disc/ |
| 9 | + |
| 10 | +# Number of attempts to make around each queued point |
| 11 | +const K := 30 |
| 12 | + |
| 13 | +## Generated points |
| 14 | +var points: PackedVector2Array |
| 15 | + |
| 16 | +var _polygon: PackedVector2Array |
| 17 | +var _bbox: Rect2 |
| 18 | + |
| 19 | +## Radius: minimum distance between generated points |
| 20 | +var _r: float |
| 21 | +## [member _r] squared – this is needed frequently so cache it for a potential |
| 22 | +## modest performance improvement. |
| 23 | +var _r_squared: float |
| 24 | + |
| 25 | +## A grid covering the area to be filled. Each element is an index into [member |
| 26 | +## points] or [code]-1[/code] if the cell is unfilled. Cells are conceptually |
| 27 | +## square, with side [member _cell_size]. |
| 28 | +var _grid: PackedInt32Array |
| 29 | + |
| 30 | +## Size of sides of cells in [member _grid], derived from [member _r]. |
| 31 | +var _cell_size: float |
| 32 | + |
| 33 | +## Points in [member _grid] are relative to this origin, so that all coordinates |
| 34 | +## are positive. |
| 35 | +var _grid_origin: Vector2 |
| 36 | + |
| 37 | +## Dimensions of [member _grid] |
| 38 | +var _grid_size: Vector2i |
| 39 | + |
| 40 | +## Elements of [member _points] that we will look to place further points near. |
| 41 | +## When we fail to place a new point near an element of this array, it is removed. |
| 42 | +## Generation is complete when this is empty. |
| 43 | +var _active: PackedVector2Array |
| 44 | + |
| 45 | + |
| 46 | +static func _bounding_box(polygon: PackedVector2Array) -> Rect2: |
| 47 | + if polygon.is_empty(): |
| 48 | + return Rect2() |
| 49 | + |
| 50 | + var bbox: Rect2 = Rect2(polygon[0], Vector2.ZERO) |
| 51 | + for point: Vector2 in polygon: |
| 52 | + bbox = bbox.expand(point) |
| 53 | + |
| 54 | + return bbox |
| 55 | + |
| 56 | + |
| 57 | +func initialise(polygon: PackedVector2Array, minimum_separation: float = 64) -> void: |
| 58 | + _polygon = polygon |
| 59 | + _bbox = _bounding_box(polygon) |
| 60 | + |
| 61 | + _r = minimum_separation |
| 62 | + _r_squared = _r * _r |
| 63 | + |
| 64 | + # Pick the cell size to be bounded by r/sqrt(2), so that each grid cell |
| 65 | + # contains at most one sample. sin(45°) = sqrt(0.5) |
| 66 | + _cell_size = _r * sqrt(0.5) |
| 67 | + _grid_origin = _bbox.position |
| 68 | + _grid_size = Vector2i( |
| 69 | + ceili(_bbox.size.x / _cell_size), |
| 70 | + ceili(_bbox.size.y / _cell_size), |
| 71 | + ) |
| 72 | + _grid.resize(_grid_size.x * _grid_size.y) |
| 73 | + assert(_grid.size() < ((1 << 31) - 1), "Grid is too large for int32 indices") |
| 74 | + _grid.fill(-1) |
| 75 | + |
| 76 | + points.clear() |
| 77 | + _active.clear() |
| 78 | + |
| 79 | + |
| 80 | +## Runs generate until no more points can be discovered |
| 81 | +func fill() -> void: |
| 82 | + while generate(): |
| 83 | + pass |
| 84 | + |
| 85 | + |
| 86 | +func _generate_first_point() -> Vector2: |
| 87 | + # TODO: Triangulate polygon; pick triangle, weighted by area; pick random point in triangle |
| 88 | + # (search keyword: barycentric). |
| 89 | + |
| 90 | + # Simpler implementation: Pick a random point on a random edge. |
| 91 | + var i := randi_range(0, _polygon.size() - 1) |
| 92 | + var j := (i + 1) % _polygon.size() |
| 93 | + return _polygon[i].lerp(_polygon[j], randf()) |
| 94 | + |
| 95 | + |
| 96 | +## Generates a single point and appends it to [member points]. |
| 97 | +## Returns [code]false[/code] when no more points can be generated. |
| 98 | +func generate() -> bool: |
| 99 | + if not points: |
| 100 | + _add_sample(_generate_first_point()) |
| 101 | + |
| 102 | + while _active: |
| 103 | + # While the active list is not empty, choose a random index |
| 104 | + # from it (say i). |
| 105 | + var n := _active.size() - 1 |
| 106 | + var i := randi_range(0, n) |
| 107 | + var p := _active[i] |
| 108 | + |
| 109 | + for j in range(K): |
| 110 | + var q := _generate_around(p) |
| 111 | + if Geometry2D.is_point_in_polygon(q, _polygon) and not _has_nearby_point(q): |
| 112 | + _add_sample(q) |
| 113 | + return true |
| 114 | + |
| 115 | + # No suitable candidate found near p; remove it from the queue |
| 116 | + _active[i] = _active[n] |
| 117 | + _active.remove_at(n) |
| 118 | + |
| 119 | + # No further points to search around |
| 120 | + return false |
| 121 | + |
| 122 | + |
| 123 | +func finished() -> bool: |
| 124 | + return _active.is_empty() |
| 125 | + |
| 126 | + |
| 127 | +func _to_grid_coords(point: Vector2) -> Vector2i: |
| 128 | + var transformed := point - _grid_origin |
| 129 | + return Vector2i( |
| 130 | + floori(transformed.x / _cell_size), |
| 131 | + floori(transformed.y / _cell_size), |
| 132 | + ) |
| 133 | + |
| 134 | + |
| 135 | +func _add_sample(point: Vector2) -> void: |
| 136 | + points.push_back(point) |
| 137 | + _active.push_back(point) |
| 138 | + |
| 139 | + var gc := _to_grid_coords(point) |
| 140 | + var grid_index := _grid_size.x * gc.y + gc.x |
| 141 | + var existing_ix := _grid[grid_index] |
| 142 | + assert( |
| 143 | + existing_ix == -1, |
| 144 | + ( |
| 145 | + "Existing point %s at grid coord %s when inserting %s" |
| 146 | + % [ |
| 147 | + points[existing_ix], |
| 148 | + gc, |
| 149 | + point, |
| 150 | + ] |
| 151 | + ) |
| 152 | + ) |
| 153 | + _grid[grid_index] = points.size() - 1 |
| 154 | + |
| 155 | + |
| 156 | +## Generates a point between r and 2r away from p, uniformly distributed. |
| 157 | +func _generate_around(p: Vector2) -> Vector2: |
| 158 | + var theta := randf_range(0, TAU) |
| 159 | + # https://stackoverflow.com/a/9048443 via https://www.jasondavies.com/poisson-disc/ |
| 160 | + var distance := sqrt(randf_range(_r_squared, 4 * _r_squared)) |
| 161 | + var q := p + (distance * Vector2.from_angle(theta)) |
| 162 | + return q |
| 163 | + |
| 164 | + |
| 165 | +## Searches [member _grid] to check if an existing point is within a distance of r from p. |
| 166 | +func _has_nearby_point(p: Vector2) -> bool: |
| 167 | + var n := 2 |
| 168 | + var gc := _to_grid_coords(p) |
| 169 | + var x0 := maxi(gc.x - n, 0) |
| 170 | + var y0 := maxi(gc.y - n, 0) |
| 171 | + var x1 := mini(gc.x + n + 1, _grid_size.x) |
| 172 | + var y1 := mini(gc.y + n + 1, _grid_size.y) |
| 173 | + for y in range(y0, y1): |
| 174 | + var gy := y * _grid_size.x |
| 175 | + for x in range(x0, x1): |
| 176 | + var existing_ix := _grid[gy + x] |
| 177 | + if existing_ix != -1: |
| 178 | + var existing := points[existing_ix] |
| 179 | + if existing.distance_squared_to(p) < _r_squared: |
| 180 | + return true |
| 181 | + return false |
0 commit comments