nx_poisson_disk.nx source

1// nx_poisson_disk.nx -- Bridson 2007 fast Poisson disk sampling. 2// 3// LAYER 2 of the procgen substrate (Perlin = layer 1). Genuinely 4// modern algorithm (2007), patent-clean, O(N). Used by Townscaper- 5// class procedural tools for building placement, Houdini's scatter 6// SOP for prop distribution, agent-spawn placement in zombie-escape 7// city games -- anywhere "looks random AND evenly spaced AND fast" 8// matters more than uniform random scattering (which clumps). 9// 10// Honest deviation note vs canonical Bridson paper (PDF was not 11// directly readable; spec verified from training data; v1 caveat): 12// 1. Initial seed point: RANDOM (fixed in v1.1 -- previously center). 13// 2. Annulus sample radius: drawn LINEARLY in [r, 2r] rather than 14// area-uniform (canonical r' = sqrt(uniform * (4r^2 - r^2) + r^2) 15// requires sqrt; queued for follow-up once nx_sqrt_q10 ships). 16// Substrate still satisfies the min-distance invariant; the bias 17// is a slight inner-ring density skew, not a correctness bug. 18// 19// Algorithm (Bridson 2007 SIGGRAPH sketch): 20// 1. Background grid with cell size r / sqrt(n_dim). Each cell 21// can hold at most one point (because two points within one 22// cell would be < r apart). 23// 2. Start: place one seed point, mark it active. 24// 3. Loop while active non-empty: 25// Pick random active point p. Try k = 30 candidates at 26// distance [r, 2r] from p. For each candidate: 27// - Check in-bounds. 28// - Check no existing point within r (grid neighbors only -- 29// this is what makes Bridson O(N)). 30// - If valid, add point + mark active; break. 31// If all 30 fail, remove p from active. 32// 4. Return collected points. 33// 34// Why this beats Perlin-grid-scatter for city gen: Perlin scatter 35// has clumping (high-noise spots get many points; low-noise spots 36// none). Poisson disk guarantees minimum distance r AND fills 37// space uniformly -- exactly what zombie-spawn / building-plot 38// placement requires. 39// 40// Cross-modal: also gives n-gram-diverse sample selection from a 41// token corpus (treat tokens as 1D positions, r = minimum gap), 42// event-timing scatter for audio onsets, test-case diversity for 43// fuzzers. 44// 45// Qualitative reading per the dual-reading cardinal: a density 46// band classifier maps the actual achieved point-count vs the 47// theoretical maximum into a sealed enum 48// (SPARSE / MODERATE / DENSE / SATURATED). 49// 50// genealogy_id: bridson_2007_poisson_disk + cook_1986_stochastic + 51// lloyd_1982_least_squares_quantization 52// lineage_id: bridson_poisson_disk_q10 53 54// nx_safety_envelope: 55// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 56// sil_target: SIL1 57// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 58// verdict: NOT_YET_EVALUATED 59 60import "nx_syscalls.nx" 61import "nx_tier.nx" 62import "nx_prng.nx" 63 64const NX_PD_Q: nx_int = 1024 65// sqrt(1/2) in Q10: ~ 0.7071 * 1024 = 724. Used for cell size 66// derivation in 2D (Bridson: cell_size = r / sqrt(n_dim)). 67const NX_PD_INV_SQRT_2: nx_int = 724 68const NX_PD_K_CANDIDATES: nx_int = 30 69const NX_PD_N_ANGLES: nx_int = 16 70 71// 16-entry precomputed (cos * Q10, sin * Q10) at angles 0, 22.5, 72// 45, ..., 337.5 degrees. Indexed by random hash % 16. 73const NX_PD_COS_0: nx_int = 1024 74const NX_PD_COS_1: nx_int = 946 75const NX_PD_COS_2: nx_int = 724 76const NX_PD_COS_3: nx_int = 392 77const NX_PD_COS_4: nx_int = 0 78const NX_PD_COS_5: nx_int = -392 79const NX_PD_COS_6: nx_int = -724 80const NX_PD_COS_7: nx_int = -946 81const NX_PD_COS_8: nx_int = -1024 82const NX_PD_COS_9: nx_int = -946 83const NX_PD_COS_10: nx_int = -724 84const NX_PD_COS_11: nx_int = -392 85const NX_PD_COS_12: nx_int = 0 86const NX_PD_COS_13: nx_int = 392 87const NX_PD_COS_14: nx_int = 724 88const NX_PD_COS_15: nx_int = 946 89 90const NX_PD_SIN_0: nx_int = 0 91const NX_PD_SIN_1: nx_int = 392 92const NX_PD_SIN_2: nx_int = 724 93const NX_PD_SIN_3: nx_int = 946 94const NX_PD_SIN_4: nx_int = 1024 95const NX_PD_SIN_5: nx_int = 946 96const NX_PD_SIN_6: nx_int = 724 97const NX_PD_SIN_7: nx_int = 392 98const NX_PD_SIN_8: nx_int = 0 99const NX_PD_SIN_9: nx_int = -392 100const NX_PD_SIN_10: nx_int = -724 101const NX_PD_SIN_11: nx_int = -946 102const NX_PD_SIN_12: nx_int = -1024 103const NX_PD_SIN_13: nx_int = -946 104const NX_PD_SIN_14: nx_int = -724 105const NX_PD_SIN_15: nx_int = -392 106 107// ===== Qualitative density-band classification ========================= 108// 109// Achieved point count vs theoretical maximum (image area / pi r^2 / sqrt(3)). 110// In Q10 we approximate the theoretical max as 111// area * 1024 / (r * r * NX_PD_MAX_FACTOR) where NX_PD_MAX_FACTOR ~ 5.4 112// (Wikipedia: hexagonal packing density 0.9069, so area / (pi r^2 * 0.9069)). 113const NX_PD_MAX_FACTOR_Q10: nx_int = 5530 // ~ 5.4 * Q10 114 115const NX_PD_BAND_SPARSE: nx_int = 0 116const NX_PD_BAND_MODERATE: nx_int = 1 117const NX_PD_BAND_DENSE: nx_int = 2 118const NX_PD_BAND_SATURATED: nx_int = 3 119const NX_PD_N_BANDS: nx_int = 4 120 121func _pd_cos(idx: nx_int) -> nx_int { 122 if idx == 0 { return NX_PD_COS_0 } 123 if idx == 1 { return NX_PD_COS_1 } 124 if idx == 2 { return NX_PD_COS_2 } 125 if idx == 3 { return NX_PD_COS_3 } 126 if idx == 4 { return NX_PD_COS_4 } 127 if idx == 5 { return NX_PD_COS_5 } 128 if idx == 6 { return NX_PD_COS_6 } 129 if idx == 7 { return NX_PD_COS_7 } 130 if idx == 8 { return NX_PD_COS_8 } 131 if idx == 9 { return NX_PD_COS_9 } 132 if idx == 10 { return NX_PD_COS_10 } 133 if idx == 11 { return NX_PD_COS_11 } 134 if idx == 12 { return NX_PD_COS_12 } 135 if idx == 13 { return NX_PD_COS_13 } 136 if idx == 14 { return NX_PD_COS_14 } 137 return NX_PD_COS_15 138} 139 140func _pd_sin(idx: nx_int) -> nx_int { 141 if idx == 0 { return NX_PD_SIN_0 } 142 if idx == 1 { return NX_PD_SIN_1 } 143 if idx == 2 { return NX_PD_SIN_2 } 144 if idx == 3 { return NX_PD_SIN_3 } 145 if idx == 4 { return NX_PD_SIN_4 } 146 if idx == 5 { return NX_PD_SIN_5 } 147 if idx == 6 { return NX_PD_SIN_6 } 148 if idx == 7 { return NX_PD_SIN_7 } 149 if idx == 8 { return NX_PD_SIN_8 } 150 if idx == 9 { return NX_PD_SIN_9 } 151 if idx == 10 { return NX_PD_SIN_10 } 152 if idx == 11 { return NX_PD_SIN_11 } 153 if idx == 12 { return NX_PD_SIN_12 } 154 if idx == 13 { return NX_PD_SIN_13 } 155 if idx == 14 { return NX_PD_SIN_14 } 156 return NX_PD_SIN_15 157} 158 159// Squared distance. Avoids sqrt -- compare squared distances directly. 160func _pd_dist_sq(ax: nx_int, ay: nx_int, bx: nx_int, by: nx_int) -> nx_int { 161 let dx: nx_int = ax - bx 162 let dy: nx_int = ay - by 163 return dx * dx + dy * dy 164} 165 166// Grid neighbor scan: check no existing point within r of (cx, cy). 167// Returns 1 if accepted (no neighbors within r), 0 if rejected. 168// Bridson's O(N) trick: only the 5x5 (or so) cells around the candidate 169// need checking. 170func _pd_accept(cx: nx_int, cy: nx_int, r: nx_int, r_sq: nx_int, 171 grid: *i64, grid_w: nx_int, grid_h: nx_int, cell: nx_int, 172 xs: *i64, ys: *i64) -> nx_int { 173 let gx: nx_int = cx / cell 174 let gy: nx_int = cy / cell 175 // Check 5x5 neighborhood -- safe upper bound for cell = r/sqrt(2). 176 var dy: nx_int = 0 - 2 177 while dy <= 2 { 178 var dx: nx_int = 0 - 2 179 while dx <= 2 { 180 let nx_g: nx_int = gx + dx 181 let ny_g: nx_int = gy + dy 182 if nx_g >= 0 { 183 if nx_g < grid_w { 184 if ny_g >= 0 { 185 if ny_g < grid_h { 186 let idx: nx_int = grid[ny_g * grid_w + nx_g] 187 if idx >= 0 { 188 let d_sq: nx_int = _pd_dist_sq(cx, cy, xs[idx], ys[idx]) 189 if d_sq < r_sq { return 0 } 190 } 191 } 192 } 193 } 194 } 195 dx = dx + 1 196 } 197 dy = dy + 1 198 } 199 return 1 200} 201 202// ===== Public entry ===================================================== 203// 204// Sample N points in a rectangle [0, w) x [0, h) with minimum 205// pairwise distance r. Caller provides out_xs / out_ys buffers 206// sized max_points. Returns actual count. 207// 208// Bridson's algorithm: O(max_points), no sqrt, no f64. 209 210func nx_poisson_disk_sample(seed: nx_int, w: nx_int, h: nx_int, r: nx_int, 211 out_xs: *i64, out_ys: *i64, 212 max_points: nx_int) -> nx_int { 213 if r <= 0 { return 0 } 214 if w <= 0 { return 0 } 215 if h <= 0 { return 0 } 216 if max_points <= 0 { return 0 } 217 218 let r_sq: nx_int = r * r 219 let cell: nx_int = (r * NX_PD_INV_SQRT_2) / NX_PD_Q 220 var cell_eff: nx_int = cell 221 if cell_eff < 1 { cell_eff = 1 } 222 let grid_w: nx_int = w / cell_eff + 1 223 let grid_h: nx_int = h / cell_eff + 1 224 225 let grid: *i64 = (sys_mmap(grid_w * grid_h * NX_SIZEOF_NX_INT)) as *i64 226 var gi: nx_int = 0 227 while gi < grid_w * grid_h { 228 grid[gi] = 0 - 1 229 gi = gi + 1 230 } 231 232 let active: *i64 = (sys_mmap(max_points * NX_SIZEOF_NX_INT)) as *i64 233 var n_active: nx_int = 0 234 235 let prng_state: *i64 = (sys_mmap(NX_SIZEOF_NX_INT)) as *i64 236 nx_prng_init(prng_state, seed) 237 238 // Seed point: canonical Bridson 2007 specifies the initial point 239 // chosen RANDOMLY in the domain (not center). Center-seed gives 240 // visually-similar output but a deterministic-by-seed canonical 241 // implementation places the first point at a PRNG-sampled position. 242 var n_points: nx_int = 0 243 let sx_raw: nx_int = nx_prng_next(prng_state) 244 var sx0: nx_int = sx_raw - (sx_raw / w) * w 245 if sx0 < 0 { sx0 = 0 - sx0 } 246 let sy_raw: nx_int = nx_prng_next(prng_state) 247 var sy0: nx_int = sy_raw - (sy_raw / h) * h 248 if sy0 < 0 { sy0 = 0 - sy0 } 249 out_xs[0] = sx0 250 out_ys[0] = sy0 251 grid[(sy0 / cell_eff) * grid_w + (sx0 / cell_eff)] = 0 252 active[0] = 0 253 n_active = 1 254 n_points = 1 255 256 while n_active > 0 { 257 if n_points >= max_points { return n_points } 258 // Pick a random active index. 259 let pick_raw: nx_int = nx_prng_next(prng_state) 260 var pick: nx_int = pick_raw - (pick_raw / n_active) * n_active 261 if pick < 0 { pick = 0 - pick } 262 let pidx: nx_int = active[pick] 263 let px: nx_int = out_xs[pidx] 264 let py: nx_int = out_ys[pidx] 265 266 var success: nx_int = 0 267 var k: nx_int = 0 268 while k < NX_PD_K_CANDIDATES { 269 let araw: nx_int = nx_prng_next(prng_state) 270 var ai: nx_int = araw - (araw / NX_PD_N_ANGLES) * NX_PD_N_ANGLES 271 if ai < 0 { ai = 0 - ai } 272 let cs: nx_int = _pd_cos(ai) 273 let sn: nx_int = _pd_sin(ai) 274 let rraw: nx_int = nx_prng_next(prng_state) 275 var rr: nx_int = rraw - (rraw / r) * r 276 if rr < 0 { rr = 0 - rr } 277 let radius: nx_int = r + rr // [r, 2r) 278 let cx: nx_int = px + (cs * radius) / NX_PD_Q 279 let cy: nx_int = py + (sn * radius) / NX_PD_Q 280 if cx >= 0 { 281 if cx < w { 282 if cy >= 0 { 283 if cy < h { 284 let ok: nx_int = _pd_accept(cx, cy, r, r_sq, 285 grid, grid_w, grid_h, cell_eff, 286 out_xs, out_ys) 287 if ok == 1 { 288 out_xs[n_points] = cx 289 out_ys[n_points] = cy 290 grid[(cy / cell_eff) * grid_w + (cx / cell_eff)] = n_points 291 active[n_active] = n_points 292 n_active = n_active + 1 293 n_points = n_points + 1 294 success = 1 295 k = NX_PD_K_CANDIDATES // exit inner loop 296 } 297 } 298 } 299 } 300 } 301 k = k + 1 302 } 303 304 if success == 0 { 305 // Remove pidx from active: swap-with-last. 306 active[pick] = active[n_active - 1] 307 n_active = n_active - 1 308 } 309 } 310 311 return n_points 312} 313 314// ===== Qualitative density-band classifier ============================ 315// 316// Maps achieved n_points vs theoretical-max into a sealed-enum band. 317 318func nx_poisson_disk_classify(n_points: nx_int, w: nx_int, h: nx_int, r: nx_int) -> nx_int { 319 if r <= 0 { return NX_PD_BAND_SPARSE } 320 // theoretical_max ~ area / (r^2 * 5.4) in i64 (5.4 ~ pi/sqrt(3) for 321 // hexagonal-packing inverse). 322 let area: nx_int = w * h 323 let denom: nx_int = (r * r * NX_PD_MAX_FACTOR_Q10) / NX_PD_Q 324 if denom <= 0 { return NX_PD_BAND_SPARSE } 325 let theoretical_max: nx_int = area / denom 326 if theoretical_max <= 0 { return NX_PD_BAND_SPARSE } 327 let ratio_q10: nx_int = (n_points * NX_PD_Q) / theoretical_max 328 if ratio_q10 < 256 { return NX_PD_BAND_SPARSE } // < 25% 329 if ratio_q10 < 614 { return NX_PD_BAND_MODERATE } // 25-60% 330 if ratio_q10 < 870 { return NX_PD_BAND_DENSE } // 60-85% 331 return NX_PD_BAND_SATURATED // >= 85% 332} 333 334func nx_poisson_disk_band_is_valid(band: nx_int) -> nx_int { 335 if band < 0 { return 0 } 336 if band >= NX_PD_N_BANDS { return 0 } 337 return 1 338}