Blue Noise as a Lattice Gibbs Ensemble

Abstract: Blue-noise sampling is widely used in computer graphics, but existing methods separate statistical modeling from scalable generation. Optimization and transport methods produce high-quality point sets by coupling all samples together. Procedural and tile-based samplers are local, but define their output only implicitly. We formulate blue-noise generation as sampling from a Gibbs distribution over binary lattice occupancies with pairwise repulsive interactions. Density, repulsion strength, interaction scale, and kernel hardness are parameters of this distribution. Because the energy sums over pairs, distant interactions can be dropped with a bounded change to the distribution, leaving a Markov random field of bounded degree. To sample it, we trace the Markov chain backward from the state we want, following Coupling Towards The Past, and cut the trace at a fixed depth. This bounds the cost, and it bounds the region each sample depends on. A tile generated on its own, with a sufficient halo, is then bit-identical to the same region generated on any larger domain, in any order and with no communication between tiles. Memory is set by the tile size, not by the output size, and accuracy is traded against cost through parameters with a proven error bound rather than by switching algorithms. We validate the model, the sampler, and these guarantees separately. The ensemble reproduces standard blue-noise spectra and moves continuously between them as its parameters vary. The sampler matches its predicted work and memory. Tiled output is verified bit-identical to full-domain generation. We demonstrate adaptive stippling at 14K, where existing methods need memory proportional to the output, along with multi-class extensions.

• The paper reformulates blue-noise generation as a binary lattice Gibbs distribution with soft pairwise repulsion, where activity, interaction strength, scale, and kernel hardness control density and spectral structure.
• The paper introduces a backward Markov-chain sampler with finite spatial range, depth, and intra-layer iterations, achieving predictable work and memory while providing a certified halo for bit-identical tiled outputs.
• Experiments show convergence by depth 5, Jacobi iteration 7, and radius 7.5, while large-scale adaptive stippling generated 1.2 million points using 50.7 MiB of peak VRAM and 975 ms per tile.

Overview

This paper reformulates blue-noise point-set generation as sampling from an explicit Gibbs distribution over binary occupancy fields on a lattice, and pairs that model with a backward Markov-chain sampler whose cost and dependency structure are bounded before any sampling begins. The work addresses a persistent split in the blue-noise literature: optimization- and transport-based methods (Lloyd relaxation, capacity-constrained Voronoi, BNOT, Gaussian blue noise) produce the best spectra but couple all samples globally, while procedural and tile-based methods execute locally but define their output only implicitly, as a procedure rather than a probability law. By defining the distribution through a pairwise energy on a lattice, the paper obtains a Markov random field of bounded degree; by tracing the chain backward from the output layer and cutting at a fixed depth, it obtains a finite dependency graph per sample. The combination yields a tiling guarantee: a tile generated with a sufficient halo is bit-identical to the same region generated on any larger domain, in any order, with no inter-tile communication. Memory is set by tile size rather than output size.

The Gibbs ensemble

The state space is $\Omega_\Lambda = \{0,1\}^{\Lambda}$ for a finite lattice domain $\Lambda \subset \mathbb{Z}^2$. The energy combines an activity term rewarding occupancy with a pairwise repulsive penalty using a generalized Gaussian kernel $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$, where $\sigma$ sets interaction scale and $p$ kernel hardness:

$$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$

The resulting distribution is $$\mu_\Lambda^\tau(v) \propto \exp[-H_{\alpha,\beta}^{\Lambda,\tau}(v)]$$, with boundary conditions handled through an exterior assignment $\tau$ chosen by the sampler. Because each close occupied pair multiplies a configuration's probability by $e^{-\beta K_{ij}} < 1$ but never forbids it, repulsion is soft: patterns remain stochastic at every finite $\beta$. At $\Lambda \subset \mathbb{Z}^2$0 the sites are independent Bernoulli variables.

The single-site conditional is available in closed form — $\Lambda \subset \mathbb{Z}^2$1, where $\Lambda \subset \mathbb{Z}^2$2 is the weighted count of occupied neighbors — which is all the sampler requires. Since the sigmoid is $\Lambda \subset \mathbb{Z}^2$3-Lipschitz, truncating the infinite-support kernel perturbs conditionals by at most $\Lambda \subset \mathbb{Z}^2$4 times the discarded tail mass, yielding a finite-range target measure with a quantified error.

The four parameters ($\Lambda \subset \mathbb{Z}^2$5, $\Lambda \subset \mathbb{Z}^2$6, $\Lambda \subset \mathbb{Z}^2$7, $\Lambda \subset \mathbb{Z}^2$8) control density, order, peak frequency, and peak sharpness respectively. Two monotonicity identities hold exactly: expected count increases with $\Lambda \subset \mathbb{Z}^2$9 (by the variance of the count), and expected pair energy decreases with $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$0 (by the variance of the pair energy). Notably, the paper argues against reading spectral order from the RAPS principal-peak height alone: raising $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$1 redistributes power across radial modes rather than amplifying the first annulus, so the scalar peak height is not an order parameter for this ensemble.

Two boundary limits connect the model to prior constructions. In the hard-core limit with a single random-order sweep, the sampler reduces exactly to dart throwing on the lattice (and Poisson-disk sampling as spacing vanishes); dart throwing is thus a boundary point of this family rather than a separate method. The other boundary, where the energy reaches its minimum, is the regime optimization methods occupy.

Spatially varying density

Adaptive generation makes both activity and scale spatially varying fields. A key technical observation is that site-dependent activities alone are insufficient: because the interaction range must match local target spacing (which scales as $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$2), fixed-scale interactions cause pairing artifacts in dense regions and under-repulsion in sparse regions. The paper therefore sets $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$3, anchored at $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$4 in the densest region and capped at $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$5, with symmetric pair weights $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$6. Activity calibration uses a local density approximation: a lookup table $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$7 precomputed on homogeneous periodic patches maps target density to activity. The LDA assumes density varies slowly relative to the local exclusion scale; rapidly varying maps incur calibration error rather than DAG failure.

Backward sampler

Rather than running a forward Gibbs chain until mixing — unbounded and hard to diagnose in strongly repulsive regimes — the sampler draws a continuous rank field ordering site updates within each sweep, plus independent coin fields, making each Gibbs update a deterministic Boolean rule. Unrolling the rank-ordered chain into a spacetime DAG indexed by depth $k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$8 gives an exact algebraic view of the chain. Three truncations make it practical:

• Spatial: retain only edges above a weight threshold ($k_{\sigma,p}(d) = \exp[-(\|d\|_2/\sigma)^p]$9), giving a finite-range target.
• Depth: keep $\sigma$0 backward layers, initializing the deepest from a product Bernoulli distribution.
• Intra-layer resolution: solve each layer's acyclic Boolean system by capped Jacobi iterations ($\sigma$1), exact if run to fixed point.

Layers are evaluated bottom-up by dynamic programming with parallel Jacobi sweeps, giving worst-case work $\sigma$2 and working memory $\sigma$3, with ranks and coins generated on the fly via counter-based RNGs. The dependency cone has deterministic radius $\sigma$4, which yields the central guarantee: with halo width $\sigma$5, tiled output is bit-identical to full-domain output for every rank field, given consistent coordinate-based seeds.

Experimental validation

Homogeneous comparisons against Lloyd, BNOT, Gaussian blue noise, dart throwing, and sliced optimal transport at matched density ($\sigma$6, $\sigma$7 periodic domain, 50 runs) place the ensemble squarely in the standard blue-noise regime. Two findings stand out. First, Lloyd achieves the highest RAPS peak but the worst anisotropy, and GBN also shows a directional peak; both persist under averaging over 50 independently initialized runs, consistent with energy minimization driving configurations toward ordered states. Ours, BNOT, and dart throwing show the flattest anisotropy curves. Second, the $\sigma$8 ablation confirms that first-peak height is not monotone in order: $\sigma$9 and $p$0 have comparable peak magnitudes under a shared vertical scale yet qualitatively different two-dimensional spectra.

Truncation ablations establish convergence thresholds: spectral statistics converge by $p$1; Jacobi solves converge by $p$2 (at $p$3 the principal peak is absent entirely); and spatial statistics converge for $p$4, capturing over 99.7% of kernel mass, with the automatic threshold retaining offsets out to $p$5. Empirical timing across a parameter grid fits $p$6 with $p$7, confirming multiplicative scaling.

The gap between certified and practical halos is substantial and honestly reported. The certified bound for the tested configuration is $p$8 pixels, but Monte Carlo tracing of 67M dependency chains finds a maximum observed radius of 39 pixels, with median 0 and $p$9 of 12.8. Disagreement between tiled and full-domain samples falls to zero empirically at $$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$0. The paper is explicit that Monte Carlo cannot rule out longer chains, so only the certified halo holds for every rank field — a factor-of-17 conservatism left uncharacterized beyond bounds.

Applications

The headline demonstration is adaptive stippling of a $$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$1 Tarantula Nebula image at 1,201,764 points, generated in $$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$2 core tiles with $$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$3 halos, bitwise identical to full-domain generation, at 50.7 MiB peak VRAM and 975 ms per tile. For comparison, GBN requires 509 MiB of VRAM (scaling linearly at ≈4.15 MiB/Mpx, $$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$4) and BNOT exceeds 16 GiB of RAM and fails at 12K resolution. This is the concrete payoff of bounded locality: constant memory regardless of output size, with tiles producible in arbitrary order. A multi-class extension replaces binary occupancy with categorical states under a Potts-style coupling $$H_{\alpha,\beta}^{\Lambda,\tau}(v) = -\alpha N_\Lambda(v) + \beta P_\Lambda^\tau(v).$$5 interpolating between per-class and combined repulsion; two-class and five-class results show well-formed blue-noise structure both per class and in the union, with the same DAG machinery unchanged.

Limitations and open questions

The paper concedes three structural costs. The lattice quantizes sample positions and confines analysis to finite domains. The pairwise energy cannot shape a spectrum the way Fattal-style density-matching energies can. Depth truncation trades exactness for predictable cost. Additionally, the sampler's applicability narrows with parameter choice: stronger order enlarges the dependency footprint (larger retained neighborhoods, more backward layers), and once that footprint approaches the working-domain size, tiling buys nothing and direct energy optimization becomes preferable. Where precisely this transition occurs is characterized only through conservative bounds; a sharper characterization would reduce the empirical nature of parameter selection. Extension to higher dimensions is identified as nontrivial, since state space, neighborhood size, and halo volume all grow quickly. Finally, the authors report that locally normalized asymmetric kernels occasionally produce visually sharper adaptive results, but such kernels admit no symmetric pairwise Hamiltonian, fall outside the analyzed model, and their stationary behavior remains an open question.

Conclusion

The paper contributes an explicit Gibbs ensemble for blue noise whose parameters continuously span a family of spectra, a backward sampler with provably bounded work, memory, and dependency cones, and a bit-exact tiling guarantee verified empirically and certified in theory. Its distinguishing property is not spectral superiority — the interior of its parameter range sits between the hard-core limit and what global optimization achieves — but the simultaneous availability of an explicit distribution, a pre-known dependency radius, and exact tiled agreement, enabling constant-memory generation at resolutions where whole-point-set methods fail outright.