# Digital Halftoning and Screening Methods for Production Printing ## From Quantization Theory to a Tileable Blue-Noise / Void-and-Cluster Threshold Array **Document class:** Engineering white paper / implementation specification **Intended audience:** developers building or maintaining halftoning/screening components, particularly in digital printing and RIP pipelines. The specification sections are written precisely enough to be handed to a coding LLM (e.g. Claude Code, Codex) as an implementation brief, but the document as a whole is meant to be read and reviewed by the engineer, not only by the tool. **Scope:** binary and multi-level halftone screening via precomputed, tileable threshold arrays, with particular focus on blue-noise spectral characteristics and the void-and-cluster construction method as the primary target algorithm. **Out of scope (mentioned only for context):** classical AM/clustered-dot screening, content-adaptive error diffusion, iterative optimization methods (Direct Binary Search) beyond a short comparative note. **© Hybrid Software Helix. Confidential document — unauthorized reproduction or distribution prohibited.** --- ## 0. How to Use This Document Sections 1–5 build the theoretical foundation: what problem we are solving, why blue noise is the right spectral target, and why void-and-cluster is a good way to construct it. A reader who only needs to *implement* the system can skip directly to Section 6 onward, but should not skip Section 5.3–5.4 (the principal-frequency relation) or Section 6.4 (why the algorithm needs three phases, not one) — these two points are the most common source of subtly-wrong implementations that "look plausible" but fail validation. Sections 6–11 are written as a direct specification: function signatures, near-Python pseudocode, parameter defaults with justification, and a concrete validation checklist with formulas. Appendix A contains a fuller reference pseudocode listing intended to be translated close to 1:1 into a real implementation language. --- ## 1. Halftoning as a Spatial Requantization Problem ### 1.1 The problem, formally Let `I(x, y)` be a continuous-tone source image with amplitude resolution of `b` bits per channel per pixel (typically `b = 8` or `16`). Let the target device support only `k` distinguishable marking states per addressable location, with `k ≪ 2^b` — in the extreme and most common case `k = 2` (mark / no mark: ink or no ink, exposed or unexposed, lit or unlit). Halftoning constructs an output image `O(x, y) ∈ {0, ..., k-1}` such that, under the spatial integration performed by the viewer's visual system (or, physically, by the finite resolving power of the eye plus any optical scattering in paper/substrate), the *local area-averaged* value of `O` over a small neighborhood approximates `I` at the corresponding location: ``` (1/|N(x,y)|) · Σ_{(x',y') ∈ N(x,y)} O(x', y') ≈ I(x, y) / (2^b - 1) ``` This is a **spatial requantization**: amplitude resolution (bits per pixel) is traded for spatial resolution (pixels per resolvable tone step). The trade is only "free" because the eye is a comparatively poor high-frequency instrument — it cannot resolve individual marks once they are small and dense enough, so it integrates them into an apparent gray level. Everything in this paper follows from taking that fact seriously as a *design constraint*, not just a convenient hand-wave. ### 1.2 Quantization noise and its spectral shaping Any per-pixel quantization to `k` levels introduces an error term: ``` e(x, y) = O(x, y) - I(x, y) / (2^b - 1) · (k - 1) ``` If `O` is produced by simple, spatially-uniform thresholding (i.e., the same fixed threshold everywhere), `e(x, y)` is strongly correlated with the image content and with itself across neighboring pixels wherever the image is locally flat — the result is **false contouring** in smooth gradients and highly visible **periodic or clustered texture** in flat midtone regions. The whole discipline of halftoning is the search for quantization strategies whose error `e(x,y)`, viewed as a 2D signal, has a *power spectrum that is inaudible to the eye* — i.e., concentrated where the visual system's sensitivity is lowest. ### 1.3 The Human Visual System as a spatial low-pass filter The relevant physiological fact is the **contrast sensitivity function (CSF)** of human vision: sensitivity to sinusoidal luminance gratings is band-pass in general, peaking at a few cycles per degree of visual angle and falling off sharply — approximately monotonically beyond the peak — at higher spatial frequencies. For the *typical* print/display viewing geometry (documents held at reading distance, or printed matter examined at arm's length), spatial frequencies above roughly 30–60 cycles/degree are effectively unresolvable; this maps to a concrete pixels-per-inch threshold once viewing distance and device resolution are fixed. Two practical consequences follow directly: - **High spatial frequencies are "free real estate."** Quantization error energy placed above the eye's cutoff is, for practical purposes, invisible. - **Low and mid spatial frequencies are expensive.** Any error energy left there — whether from a fixed threshold, an unlucky periodic dither matrix, or a poorly-diffused error-diffusion kernel — will be visible as texture, banding, or moiré. This is precisely the design target formalized as "blue noise" in Section 5: push all of the quantization error's spectral energy above a frequency the HVS cannot resolve, and make sure none of it concentrates into narrow, high-energy spectral peaks (which the eye is *very* good at detecting even when they're nominally "high frequency," because periodicity itself is a strong perceptual cue independent of the CSF envelope). ### 1.4 Two structurally different strategies 1. **Spatially-varying, precomputed threshold (ordered dithering).** Compare each pixel to a threshold that varies by position according to a fixed, precomputed matrix — a *dither array* / *threshold array* / *screen*. The matrix can be periodic (classical AM/clustered-dot screening, Bayer dithering) or aperiodic-looking (stochastic/blue-noise screening, the subject of this paper). 2. **Feedback quantization (error diffusion).** Use a trivial fixed threshold, but *diffuse* the resulting quantization error into not-yet-processed neighboring pixels, so that the running local average tracks the source image. Content-adaptive, sequential, not precomputable as a static tile. This paper's deliverable belongs to family (1): a **precomputed, tileable, aperiodic threshold array with blue-noise spectral properties**, generated via the **void-and-cluster** method. Family (2) is covered in Section 4 only for context and contrast, because understanding *why* error diffusion is not the right tool for a tile-based production pipeline is useful background for anyone maintaining this system. --- ## 2. A Brief, Contextual Note on AM (Clustered-Dot) Screening Classical **amplitude-modulated (AM) screening** is the photomechanical-era solution and remains the default in much of offset and commercial print production. Tone is encoded by dot *size* on a periodic grid at a fixed **screen ruling** (lines per inch, lpi) and **screen angle**. For color reproduction, each ink separation is screened at a different angle — classically 15°/75°/0°/45° for C/M/Y/K — specifically to control the interference pattern (**moiré**) that inevitably arises when multiple periodic gratings of similar frequency are superimposed; the angle offsets push the resulting beat pattern (the **rosette**) into a frequency and orientation that is comparatively unobjectionable, rather than eliminating it. Two structural properties of AM screening matter for contrast with this paper's target method: - It is **periodic by construction** — its power spectrum is a discrete comb of harmonics at the screen frequency and its multiples, not a continuous band. This is exactly what makes moiré, and interference with the marking engine's native addressability or with scanning/imaging pitches, a first-order concern that must be engineered around (angle selection, rational tangent screening, supercell techniques). - It modulates **cluster size**, not dot position: at a given screen cell, tone is expressed as *how large* the single dot cluster is, not *where* dots are placed. AM screening is not treated further in this document. It is mentioned so the implementer understands what problem stochastic/blue-noise screening is *not* trying to replicate: this paper's method is aperiodic (or periodic only at a tile size chosen to be visually and mechanically irrelevant) and modulates dot **position/density** at close to the addressable-pixel scale, which is what makes cross-channel moiré manageable through decorrelation (Section 9) rather than through angle geometry. --- ## 3. Ordered Dithering: Threshold Arrays and Dither Matrices ### 3.1 Formal definition Ordered dithering compares each source pixel value against a threshold drawn from a small `W × H` matrix, tiled across the image with wraparound indexing: ``` O(x, y) = ON if I(x, y) > T( x mod W, y mod H ) = OFF otherwise ``` Following Ulichney's terminology, it is useful to separate two related objects: - A **dither pattern** stores *threshold indices*: an ordering, `0..(W·H − 1)`, of the tile's positions by "priority" — the sequence in which positions turn ON as the (uniform, hypothetical) input tone rises from 0 to maximum. This is a pure combinatorial object, independent of bit depth. - A **threshold array** is the dither pattern rescaled to the tonal range of the actual source image (Section 5.1's rank-to-threshold formula), and is the object actually compared against pixel values at runtime. This separation matters in practice: the *rank matrix* (dither pattern) is what the void-and-cluster algorithm computes and what should be stored/versioned as the canonical artifact; the threshold array is a cheap, reproducible rescaling of it, and — as argued in Section 10 — should stay separate from any device-specific tone curve. ### 3.2 The Bayer matrix, and why it is not the target here The classical **Bayer matrix** is constructed recursively: ``` M_1 = [0] M_{2n} = [ 4·M_n 4·M_n + 2 ] [ 4·M_n + 3 4·M_n + 1 ] ``` producing, for `n=2`, the familiar 4×4 array. Its recursive, self-similar structure makes it trivial to generate and gives it a hard mathematical guarantee of a well-defined threshold at every scale — but that same recursive structure means its power spectrum is concentrated in a handful of discrete low-order harmonics of the tile frequency, not spread as a continuous high-frequency band. The visible symptom is the well-known cross-hatch/checkerboard texture of Bayer-dithered midtones, and — more importantly for this document's purposes — badly *uneven* quality across the tonal range: some gray levels look tolerable, others (particularly near 25%/75%, and at very light/dark tones) show strong, obviously periodic structure. Bayer dithering is legitimate as a cheap, fully deterministic fallback for extremely constrained real-time contexts (embedded displays with no room for a large precomputed LUT), but it is explicitly **not** the target of this specification. It is mentioned so the implementer understands the failure mode that motivates everything that follows: a "correct-looking" dither matrix can still be spectrally bad in a way that's easy to miss by eye at a single tone level and only becomes obvious across the full ramp, or under the radially-averaged power spectrum test of Section 11. ### 3.3 What makes a dither array "good" Ulichney's criteria — which this paper adopts as the acceptance bar for the deliverable — for a well-formed dither array are: 1. **Homogeneity at every gray level**, not just at 50%: the set of ON pixels below any threshold should be, itself, a well-distributed pattern. 2. **No low-frequency energy** below a level-dependent principal frequency (Section 5.2). 3. **A flat, isotropic (direction-independent) high-frequency spectrum** — genuine noise-like structure, free of both periodic peaks and directional bias. 4. **Seamless tileability** — no visible discontinuity where one tile repetition meets the next. Requirement 1 is exactly why the void-and-cluster algorithm is built as a monotonic full-tile ranking (Section 6) rather than as an optimization of a single tone level. Requirements 2–3 are the formal definition of blue noise (Section 5). Requirement 4 is an engineering constraint specific to *this document's* deployment target (a reusable production tile) and is what drives the toroidal-convolution requirement in Section 6.3. --- ## 4. Error Diffusion: Context and Contrast Error diffusion is covered here only to give the implementer a clear mental model of the alternative family, and to justify why it is excluded from this specification's deliverable. ### 4.1 Mechanism For each pixel, in scan order: 1. Threshold the pixel (using the *original* source value plus any error already diffused into it from prior pixels) against a single fixed threshold (typically the midpoint). 2. Compute the quantization error: the difference between the thresholded output and the (adjusted) input value. 3. Distribute that error, weighted by a fixed kernel, into neighboring *not-yet-processed* pixels. ### 4.2 Classical kernels | Method | Kernel (weights, denominator shown) | Notes | |---|---|---| | Floyd–Steinberg (1976) | `· * 7 / 5 3 1` (÷16) | The original; simple, fast, prone to "worm" artifacts in some tone ranges. | | Jarvis, Judice & Ninke (1976) | 12-tap, 2-row-ahead kernel (÷48) | Smoother, more diffuse error spread; slower. | | Stucki (1979) | 12-tap variant of JJN (÷42) | Sharper than JJN, fewer artifacts than Floyd–Steinberg. | | Sierra / Sierra Lite | compact 2–3 tap variants | Cheap approximations for constrained environments. | | Ostromoukhov (2001) | tone-dependent variable coefficients, looked up per input level | State of the art among error-diffusion kernels; matches blue-noise quality across nearly the full tonal range without the specific artifacts of fixed kernels. | Serpentine (boustrophedon) scan order — alternating left-to-right and right-to-left on successive rows — is a standard refinement that reduces the directional bias that a fixed raster scan order otherwise introduces. ### 4.3 Why error diffusion is out of scope here Error diffusion is genuinely excellent at per-pixel image quality, and (with tone-dependent kernels) can itself achieve blue-noise-like spectral characteristics as an emergent property of its feedback structure. But it has three properties that make it the wrong tool for *this specification's* deliverable: - **Content-dependence.** The output pattern depends on the actual image content, not just on position — so it cannot be precomputed once and reused as a static tile/LUT the way a threshold array can. - **Sequential data dependency.** Each pixel's decision depends on the accumulated error from prior pixels in scan order, which resists straightforward parallelization (GPU-friendly variants exist — block-based, or the "line diffusion" family — but they are nontrivial adaptations, not the base algorithm). - **Directional artifacts.** Even with serpentine scanning, some tone ranges show visible low-frequency directional texture ("worms") that a well-constructed blue-noise threshold array does not exhibit. If a future requirement calls for content-adaptive, non-tile-based screening (e.g. a high-quality soft-proofing path where GPU throughput is not the binding constraint), error diffusion — specifically the Ostromoukhov variable-coefficient variant — is the natural next thing to implement, and can reuse the same validation methodology from Section 11. It is explicitly not what Section 6 specifies. --- ## 5. Blue Noise: Formal Definition and Spectral Requirements ### 5.1 The spectral definition A binary dot pattern `P` (at a given fraction of ON pixels, `g`) is said to exhibit **blue-noise characteristics** (Ulichney, 1987/1988) if its **radially-averaged power spectrum (RAPS)** — the 2D power spectrum of `P`, averaged over all directions at each spatial-frequency radius — satisfies: 1. **No energy at low frequency**, below a level-dependent **principal frequency** `f_g` (defined in 5.2). This encodes the *absence of large-scale structure*: no visible clumping, no striping, no macroscopic gradient artifacts. 2. **A sharp transition** at `f_g`, from "no energy" to "full energy" — not a gradual ramp. 3. **A flat, isotropic spectrum above `f_g`.** "Flat" (in radius) means no dominant spatial frequency — genuine noise, not periodicity (which shows as a sharp peak in the RAPS). "Isotropic" (no directional bias) means the 2D spectrum, not just its radial average, is close to rotationally symmetric — a pattern can have a radially "flat-looking" RAPS while still being strongly anisotropic (e.g. concentrated on the diagonals), which would still be visible as directional texture; the isotropy check is therefore a *separate* diagnostic from the RAPS radial shape (see Section 11.3). ### 5.2 The principal frequency `f_g` is derived from the statistics of a maximally-disordered (Poisson-disk-like) point process at density `g`: if `N` points are placed with no two closer than the "natural" minimum spacing implied by their density, the expected nearest-neighbor spacing scales as `1/sqrt(g)` (in pixels, for `g` = fraction ON per unit area), giving: ``` f_g ≈ sqrt(g) for g ≤ 0.5 f_g ≈ sqrt(1 - g) for g > 0.5 (by symmetry: treat OFF pixels as the "dots" above 50%) ``` This is the frequency at which the dither pattern's structure should be concentrated at density `g`: not zero (that would imply either total disorder down to arbitrarily fine scale, i.e. white noise, or systematic clustering) and not the Nyquist limit (that would only be achievable, and only meaningfully define "blue," at `g = 0.5`). ### 5.3 Blue noise vs. neighboring "noise colors" The color-of-noise vocabulary, borrowed from the light-spectrum analogy, is used loosely but consistently in this literature: - **White noise**: flat spectrum at *all* frequencies including low ones — i.e., unconstrained random dot placement. This is what you get from naive per-pixel random thresholding, and it looks grainy/clumpy because, by chance, it contains both large voids and tight clusters. - **Blue noise**: as defined above — no low-frequency energy, flat/isotropic above `f_g`. This is the target. - **Green noise**: mid-frequency-concentrated, associated with dot clustering at a controlled coarser scale — relevant to stochastic *clustered-dot* screening (a hybrid between AM and blue-noise approaches) but not this paper's target. - **Violet/purple noise**: energy increasing with frequency (as opposed to blue noise's roughly flat high-frequency plateau) — occasionally discussed in the error-diffusion literature; not directly relevant here. ### 5.4 Why anisotropy is a distinct failure mode A dither array can pass a naive "radially averaged spectrum looks flat above `f_g`" check while still being visually poor, if the *directional* distribution of that energy is uneven — e.g. concentrated along the horizontal/vertical axes (a common artifact of poorly-tuned or non-toroidal implementations) or along the diagonals. This shows up perceptually as faint grid- or diagonal-hatch texture even when no single frequency dominates radially. Section 11.3 gives a concrete anisotropy diagnostic; the implementer should not treat "RAPS looks flat" as sufficient evidence of quality on its own. ### 5.5 The Fourier-domain alternative: Blue Noise Masks Mitsa & Parker (1991/1992) construct dither arrays by iteratively filtering each level's dot profile in the *frequency domain* with an explicit blue-noise-shaped filter, then thresholding back to a binary pattern and proceeding level-by-level (adding/removing minimal numbers of dots between adjacent levels to preserve the monotonic-ranking property). This produces very high-quality masks and gives more explicit control over the target spectral shape (since the filter is specified directly), but it is computationally heavier and materially more complex to implement correctly (repeated FFTs per level, careful handling of the level-to-level dot-count transition) than the purely spatial-domain void-and-cluster method. **This specification does not target the Blue Noise Mask approach**; it is noted here so the implementer understands it as the higher-cost, higher-explicit-control alternative to void-and-cluster, should quality requirements ever exceed what void-and-cluster delivers. --- ## 6. Void-and-Cluster: Concept and Full Derivation ### 6.1 What the algorithm actually optimizes Void-and-cluster (Ulichney, 1993) builds a full **rank matrix** — a permutation of every position in an `N × N` tile, assigning each position an integer rank `0..N²−1` — such that **every prefix** of that ranking (i.e., the set of positions with rank below any cutoff, for every possible cutoff) is, on its own, a well-distributed blue-noise binary pattern at that density. This is strictly stronger than optimizing a single density (e.g. 50%): it guarantees the monotonic, non-flickering, non-clumping progression from black to white needed for a reusable production screen, because in production the same tile is thresholded at every gray level present in the actual image, not just one. ### 6.2 The two primitive operations Both primitives are defined via a **Gaussian-filtered energy map** `E = P ⊛ G_σ`, where `P` is the current binary pattern (1 for ON, 0 for OFF) and `G_σ` is a 2D isotropic Gaussian kernel with standard deviation `σ`, and `⊛` denotes **toroidal (wrap-around) convolution** — see 6.3 for why this is non-negotiable. - **Tightest cluster**: `argmax_{(x,y): P(x,y)=1} E(x,y)` — among currently-ON pixels, the one sitting in the locally densest neighborhood of other ON pixels. Removing this pixel does the least damage to the pattern's overall distribution (it was the most "redundant" ON pixel). - **Largest void**: `argmin_{(x,y): P(x,y)=0} E(x,y)` — among currently-OFF pixels, the one farthest from any ON pixel. Adding a pixel here does the most to homogenize the pattern (it fills the biggest gap). Intuitively: `E` is a smoothed density map of "how ON-ish is the neighborhood around this point." High `E` at an ON pixel = redundant cluster member. Low `E` at an OFF pixel = center of an empty region. ### 6.3 The Gaussian kernel and toroidal wrap **Kernel shape.** A 2D isotropic Gaussian is used because it is smooth, rotationally symmetric (so it does not itself introduce directional bias — see 5.4), and has a single, easily-tuned scale parameter `σ`. `σ ≈ 1.5` pixels is Ulichney's original recommendation and is a reasonable default; the practical tuning range is roughly `1.0–2.5` (see Section 8.1 for the trade-off). **Truncation.** A Gaussian's contribution falls below numerical/visual significance beyond about `3σ`; in practice, precompute the kernel as a finite `⌈6σ⌉ × ⌈6σ⌉` (or `⌈6σ⌉+1`, to keep an odd, centered window) patch rather than convolving over the whole tile. **Toroidal wrap — mandatory, not optional.** Every distance and every convolution must wrap around the tile edges: for two coordinates `a, b` on an axis of size `N`, the toroidal distance is `min(|a-b|, N-|a-b|)`, and the Gaussian weight is computed from that wrapped distance on each axis independently before combining (or equivalently, the convolution is a circular convolution, `mod N` on both axes). This is precisely what guarantees the final tile abuts itself with **no visible seam** when repeated across a full page or plate. Skipping this — e.g. using a clamped-at-the-edge or zero-padded convolution, which is what most off-the-shelf image-filtering library defaults will silently do — is the single most common implementation defect: the resulting mask can look fine evaluated on its own tile and only reveal a visible discontinuity once actually tiled at production scale. **Any implementation that calls a generic library convolution function must explicitly verify it is using circular/wrap boundary conditions, not the library's default.** ### 6.4 The three phases, precisely **Phase 0 — Initial balanced pattern (the "prototype").** Generate a random binary pattern with a small initial count `n0` of ON pixels (roughly 10% of the tile area is Ulichney's suggested starting point; the exact fraction is not critical to final quality, only to Phase-0 convergence time). Then iterate: ``` loop: c = tightest_cluster(P) # among ON pixels v = largest_void(P) # among OFF pixels if c == v: # same position selected both ways → converged break P[c] = OFF P[v] = ON ``` This converges to a locally-optimal, homogeneously-distributed configuration for exactly `n0` ON pixels — the **prototype binary pattern** — which seeds both ranking phases below. (Convergence is empirical, not formally guaranteed to terminate in bounded steps for all inputs, but is reliable in practice for reasonable `σ` and `n0`; cap the loop with a generous iteration limit — e.g. `10 × n0` — as a safety net, and treat hitting the cap as a signal to inspect parameters rather than as a hard failure.) **Phase 1 — Rank the initial minority (ranks `0 .. n0−1`), counting down.** Starting again from the prototype (with `n0` ON pixels), repeatedly: ``` remaining = copy(prototype) for rank in range(n0 - 1, -1, -1): c = tightest_cluster(remaining) rank_matrix[c] = rank remaining[c] = OFF ``` This produces the correct ordering for the **darkest** end of the tonal range: the very first pixels to switch ON as gray level rises from 0 are exactly the ones identified here as "least essential to the prototype cluster" — i.e., the ones that, when the prototype was built up from nothing, would have been the last, most marginal additions. Removing them one at a time from most-to-least redundant and assigning descending ranks reconstructs that order directly. **Phase 2 — Rank from `n0` up to `N²/2 − 1`, counting up, by void-filling.** Starting again from the prototype: ``` building = copy(prototype) for rank in range(n0, N*N // 2): v = largest_void(building) rank_matrix[v] = rank building[v] = ON ``` Each step finds the currently-largest gap and fills it, which is exactly what keeps the pattern homogeneous as density increases past the prototype's own density, up to the halfway point. **Phase 3 — Rank from `N²/2` up to `N² − 1`, counting up, by inverse-cluster-filling.** For the upper half of the tonal range the roles invert: work on the **complement** pattern (swap ON/OFF) and repeatedly find the tightest cluster **of OFF pixels in the complement view** — equivalently, the position that, if turned ON in the real pattern, would most effectively fill in the tightest remaining cluster of still-empty positions: ``` complement = invert(building) # 'building' now has N²/2 ON pixels from Phase 2 for rank in range(N*N // 2, N*N): c = tightest_cluster(complement) # tightest cluster of OFF-in-real / ON-in-complement... # concretely: find the OFF position in 'building' most surrounded by OTHER OFF positions rank_matrix[c] = rank building[c] = ON complement[c] = OFF ``` (Implementation note: Phase 3's "tightest cluster of the OFF set" is computed with the same `tightest_cluster` primitive as Phase 1, applied to the *complement* pattern — no new primitive is needed, only a re-use of the existing one against inverted data.) ### 6.5 Why three phases, and not one pass in each direction A tempting simplification is to run only a single void-filling pass from empty to full (skip Phase 1, and treat Phase 3 as "more of Phase 2"), or symmetrically, only a single cluster-removal pass from full to empty. Both simplifications fail the homogeneity-at-every-level requirement (Section 3.3, requirement 1) in a specific, predictable way: - **"Best void to fill next"** and **"least essential existing cluster member"** are genuinely different criteria, and they diverge specifically near the tonal extremes (very light and very dark), where the number of ON (or OFF) pixels is small and the two criteria's rankings of the "next" pixel to touch can disagree. - Phase 1's role is specifically to get the *darkest* tones right by working backward from a pattern that is already known-good at density `n0`, rather than forward from nothing (which would re-derive an ordering with no such known-good anchor point). - Phase 3's role is the mirror-image problem at the *lightest* tones: by density `N²/2`, continuing to "find the largest void" among an ever-shrinking OFF set becomes numerically and qualitatively a worse-conditioned problem than switching to "find the tightest cluster of the complement," which is the same well-conditioned problem as Phase 1, just applied to the inverted pattern. Empirically (per Ulichney 1993 and the substantial literature building on it), this three-phase construction is what delivers blue-noise-quality RAPS characteristics across the *entire* 0–100% range, not only near the density used to build the prototype. A single-pass simplification will typically look acceptable at mid-tones (where it's closest to the prototype's own density) and visibly worse — often reverting toward white-noise-like graininess or unwanted micro-clustering — near the shadows or highlights, which is exactly where a production RIP most needs a clean screen (shadow/highlight dot control is a classic hard problem in physical print reproduction, independent of the screening algorithm). ### 6.6 Relationship to Direct Binary Search (DBS) **Direct Binary Search** (Allebach and collaborators) is a higher-quality, substantially more expensive alternative: instead of the greedy, single-pass construction above, DBS iteratively perturbs a binary pattern (swap, toggle) and accepts each change only if it improves an explicit human-visual-system-weighted error metric against the ideal continuous-tone target, iterating to convergence (often with simulated-annealing-style acceptance). DBS-generated screens are generally considered the highest-quality achievable and are the subject of most recent literature extensions (clustered-dot DBS variants, multiscale/multipass refinements, and — most recently — learned/reinforcement-learning approaches that attempt to approximate DBS quality at a fraction of the computational cost). **DBS is out of scope for this specification** — it is orders of magnitude more expensive to compute (iterative optimization over the full pattern, evaluated against an HVS model, versus void-and-cluster's single greedy pass) — but it is the natural next step if a future quality bar exceeds what void-and-cluster delivers, and it is worth the implementer knowing the name for that reason. --- ## 7. Implementation Brief ### 7.1 Deliverables ``` generate_prototype(size, seed_fraction=0.10, sigma=1.5, rng_seed=None) -> BinaryPattern void_and_cluster_rank(size, sigma=1.5, seed_fraction=0.10, rng_seed=None) -> RankMatrix rank_to_threshold(rank_matrix, bit_depth=8) -> ThresholdArray apply_screen(image, threshold_tile) -> BinaryImage validate_blue_noise(threshold_tile, levels=[...]) -> ValidationReport ``` `size` must be a positive integer (square tiles recommended; see 8.2 for non-square considerations). Prefer powers of two for `size` (Section 8.2). ### 7.2 The incremental energy map (critical for performance) The naive implementation recomputes the full convolution `P ⊛ G_σ` from scratch every time a pixel is toggled. That is correct but wasteful: across Phase 0's convergence loop plus Phases 1–3, the pattern is toggled roughly `O(N²)` times (plus Phase 0's iterations), and a full `N×N` convolution per toggle makes the whole generation process scale far worse than necessary. **Maintain a persistent `energy[N][N]` array and update it incrementally.** Toggling a single pixel changes `E` only in a local neighborhood: precompute the truncated Gaussian kernel patch once (a small `(6σ+1) × (6σ+1)` array), and on each toggle, add (if turning ON) or subtract (if turning OFF) that kernel patch into `energy`, centered at the toggled position, using toroidal (wrapped) indexing for the patch's footprint. This reduces each toggle's cost from `O(N²)` to `O(σ²)`, which is what makes generating tiles up to several hundred pixels on a side computationally practical. ```python # Reference-quality pseudocode. Translate directly; do not "simplify away" the toroidal wrap. import numpy as np def make_gaussian_patch(sigma): radius = int(np.ceil(3 * sigma)) ax = np.arange(-radius, radius + 1) xx, yy = np.meshgrid(ax, ax) kernel = np.exp(-(xx**2 + yy**2) / (2 * sigma**2)) return kernel, radius def toggle_pixel(energy, kernel, radius, x, y, N, turning_on: bool): sign = 1.0 if turning_on else -1.0 for dy in range(-radius, radius + 1): for dx in range(-radius, radius + 1): wx = (x + dx) % N wy = (y + dy) % N energy[wy, wx] += sign * kernel[dy + radius, dx + radius] def tightest_cluster(pattern, energy): # argmax of energy restricted to ON pixels masked = np.where(pattern == 1, energy, -np.inf) idx = np.unravel_index(np.argmax(masked), masked.shape) return idx # (y, x) def largest_void(pattern, energy): # argmin of energy restricted to OFF pixels masked = np.where(pattern == 0, energy, np.inf) idx = np.unravel_index(np.argmin(masked), masked.shape) return idx # (y, x) ``` Finding the current tightest cluster / largest void is a masked argmax/argmin over `energy` — an `O(N²)` scan per rank assigned in the straightforward implementation shown above. For the tile sizes recommended in this document (up to a few hundred pixels per side), a plain masked-scan per step is fast enough to run as an offline precomputation step; do not add a maintained priority-queue/heap structure unless profiling on the actual target tile size shows the scan is the bottleneck — it adds real implementation complexity (keeping a heap consistent under the incremental energy updates) for a cost that is, at these tile sizes, usually not the dominant term. ### 7.3 Full generation pseudocode ```python def void_and_cluster_rank(size, sigma=1.5, seed_fraction=0.10, rng_seed=None): rng = np.random.default_rng(rng_seed) N = size kernel, radius = make_gaussian_patch(sigma) # --- Phase 0: build the prototype --- n0 = max(1, int(round(seed_fraction * N * N))) pattern = np.zeros((N, N), dtype=np.uint8) flat_idx = rng.choice(N * N, size=n0, replace=False) pattern.flat[flat_idx] = 1 energy = np.zeros((N, N), dtype=np.float64) for idx in flat_idx: y, x = divmod(int(idx), N) toggle_pixel(energy, kernel, radius, x, y, N, turning_on=True) max_iters = 10 * n0 for _ in range(max_iters): cy, cx = tightest_cluster(pattern, energy) vy, vx = largest_void(pattern, energy) if (cy, cx) == (vy, vx): break pattern[cy, cx] = 0 toggle_pixel(energy, kernel, radius, cx, cy, N, turning_on=False) pattern[vy, vx] = 1 toggle_pixel(energy, kernel, radius, vx, vy, N, turning_on=True) prototype = pattern.copy() proto_energy = energy.copy() rank_matrix = np.full((N, N), -1, dtype=np.int64) # --- Phase 1: rank 0 .. n0-1, counting down, removing tightest clusters --- work = prototype.copy() work_energy = proto_energy.copy() for rank in range(n0 - 1, -1, -1): cy, cx = tightest_cluster(work, work_energy) rank_matrix[cy, cx] = rank work[cy, cx] = 0 toggle_pixel(work_energy, kernel, radius, cx, cy, N, turning_on=False) # --- Phase 2: rank n0 .. N*N/2 - 1, counting up, filling largest voids --- work = prototype.copy() work_energy = proto_energy.copy() half = (N * N) // 2 for rank in range(n0, half): vy, vx = largest_void(work, work_energy) rank_matrix[vy, vx] = rank work[vy, vx] = 1 toggle_pixel(work_energy, kernel, radius, vx, vy, N, turning_on=True) # --- Phase 3: rank N*N/2 .. N*N - 1, counting up, filling tightest OFF-clusters --- complement = 1 - work # OFF-in-'work' becomes ON-in-complement comp_energy = np.zeros((N, N), dtype=np.float64) # rebuild energy for complement fresh (cheap: O(N^2), done once) on_positions = np.argwhere(complement == 1) for (yy, xx) in on_positions: toggle_pixel(comp_energy, kernel, radius, int(xx), int(yy), N, turning_on=True) for rank in range(half, N * N): cy, cx = tightest_cluster(complement, comp_energy) rank_matrix[cy, cx] = rank complement[cy, cx] = 0 toggle_pixel(comp_energy, kernel, radius, cx, cy, N, turning_on=False) assert np.all(rank_matrix >= 0), "every position must be ranked exactly once" assert sorted(rank_matrix.flatten().tolist()) == list(range(N * N)), \ "rank matrix must be a valid permutation of 0..N^2-1" return rank_matrix def rank_to_threshold(rank_matrix, bit_depth=8): N2 = rank_matrix.size max_level = (1 << bit_depth) - 1 normalized = (rank_matrix.astype(np.float64) + 0.5) / N2 return np.round(normalized * max_level).astype( np.uint8 if bit_depth <= 8 else np.uint16 ) def apply_screen(image, threshold_tile): H, W = image.shape th, tw = threshold_tile.shape yy, xx = np.indices((H, W)) tiled_threshold = threshold_tile[yy % th, xx % tw] return (image > tiled_threshold).astype(np.uint8) ``` Rebuilding the complement's energy map from scratch at the Phase 2→3 boundary (rather than incrementally deriving it from `work_energy`) is deliberate: it is a single `O(N²)`-cost operation performed exactly once per generation run, it sidesteps any risk of accumulated floating-point drift across ~`N²` incremental updates, and it keeps the Phase 3 code path independent and easy to test in isolation. ### 7.4 Complexity summary - Kernel precompute: `O(σ²)`, once. - Phase 0: up to `max_iters` toggles, each `O(σ²)` for the energy update and `O(N²)` for the masked argmax/argmin scan → `O(max_iters · N²)` dominated by the scan. - Phases 1–3 together: `N²` rank assignments, each `O(σ²) + O(N²)` → `O(N⁴)` dominated by the scan, plus one `O(N²)` complement-energy rebuild. - **Net: `O(N⁴)`** for the reference (masked-scan) implementation. For `N = 256` this is on the order of `4×10⁹` elementary comparisons — seconds to low tens of seconds in a reasonably optimized native/vectorized implementation, and this is an **offline, one-time precomputation** per tile configuration, not a runtime cost. If a future requirement pushes `N` well beyond 256–512 and precompute time becomes inconvenient, the standard optimization is replacing the masked linear scan with a maintained max-heap / min-heap (or a coarser spatial index) over `energy`, bringing the per-step search down toward `O(log N²)` at the cost of nontrivial bookkeeping to keep the heap consistent under the incremental energy updates — implement this only if profiling on the real target size demonstrates the need. ### 7.5 Determinism and reproducibility Always accept and propagate an explicit `rng_seed`. Given a fixed `(size, sigma, seed_fraction, rng_seed)`, generation must be **bit-for-bit reproducible** — this matters for regression testing (Section 11) and for reproducing a specific production tile from a stored parameter set rather than having to archive the (potentially large) rank matrix itself, though archiving the final rank matrix as the canonical shipped artifact is still recommended (Section 10.1) since regeneration cost is nontrivial (Section 7.4). --- ## 8. Parameter Tuning Guide ### 8.1 Gaussian `σ` - **Too small** (`σ ≲ 0.8`): the filter under-corrects longer-range clustering; the resulting pattern can retain visible low-frequency structure, failing the "no energy below `f_g`" requirement. - **Too large** (`σ ≳ 3`): computation cost per toggle grows as `σ²`; beyond a point, quality does not improve and can mildly degrade, because the algorithm starts optimizing against correlations at a scale larger than what actually matters for the target print/display resolution. - **Default: `σ = 1.5`** (Ulichney's original value), tunable in the roughly `1.0–2.5` range. If a specific target resolution and viewing distance are known, it is reasonable to tie `σ` loosely to the expected dot pitch, but empirical validation (Section 11) should be the final arbiter, not a formula. ### 8.2 Tile size `N` - **Larger `N`** → longer effective spatial period before the pattern repeats → less perceptible large-scale repetition, but linearly more memory for the stored LUT and — per Section 7.4 — roughly `O(N⁴)` more precompute time. - **Recommended production range: 64×64 to 256×256.** 256×256 is a solid default for high-quality output (large-format print, fine inkjet); 128×128 is a good quality/cost compromise; 64×64 is adequate for screen preview, fast proofing, or lower-resolution engines. - **Prefer powers of two.** This is not required by the algorithm itself, but it makes the runtime tiling operation in `apply_screen` a cheap bitmask (`x & (N-1)` instead of `x % N`), which matters when this function sits in a real-time RIP/rasterization hot path (relevant given a GPU-accelerated RIP pipeline as the eventual deployment target). - **Non-square tiles** are possible (independent `sigma_x`/`N_x` and `sigma_y`/`N_y`) but add complexity for a benefit that is rarely worth it in practice; default to square tiles unless a specific device aspect-ratio constraint demands otherwise. ### 8.3 Seed fraction (Phase 0) Controls only the starting point and convergence speed of Phase 0, not the final quality of the rank matrix in any first-order way (the whole point of Phase 0's iterate-to-convergence loop is to erase dependence on the initial random draw). `0.10` (10% of tile area) is a reasonable, well-tested default; values from `0.05` to `0.20` behave similarly. Extremely small (`<1%`) or large (`>40%`) values are not recommended: they push Phase 0 away from the density range where its convergence behavior is well characterized in the literature. ### 8.4 Bit depth of the threshold array Match the precision needs of the pipeline: 8-bit (`0–255`) is sufficient for most direct binary/print-engine use; use 16-bit if the threshold array will be composed with a high-precision tone-response curve (Section 10) before final quantization, to avoid introducing additional banding from the threshold array's own rounding. --- ## 9. Multi-Level and Multi-Channel Extension ### 9.1 Multi-level (multi-drop) output If the target device supports more than two states per pixel (e.g., inkjet with light/medium/heavy drop sizes, or a multi-bit exposure engine), the same rank matrix generalizes directly: instead of one threshold, define `k − 1` cut points along the `0..N²−1` rank range (evenly spaced, or perceptually/device-spaced if the device's own level-to-optical-density mapping is non-uniform) to produce `k` output levels. Because the rank matrix already encodes a single monotonic, blue-noise-consistent ordering across the *entire* density range by construction (Section 6.4's whole point), each of the resulting `k − 1` binary transitions individually inherits the blue-noise property — no separate mask needs to be generated per level. ### 9.2 Multi-channel decorrelation For independent color channels (CMYK or other multi-ink systems), naively reusing the *same* rank matrix for every channel will produce strongly correlated dot patterns across channels — visually, this shows up as coincident dot placement (worse graininess than either channel alone) at some tone combinations and visible color moiré at others, because a shared spatial pattern, thresholded at different levels for different channels, does not decorrelate simply by virtue of the levels differing. Generate **independently-seeded masks per channel** (different `rng_seed` per channel is the simplest approach and is generally sufficient) so that each channel's dot pattern is blue-noise-correct on its own but spatially uncorrelated with the others. This is the direct, stochastic-screening analogue of AM screening's angle separation (Section 2) — achieved through decorrelated randomness in the mask-generation seed rather than through geometric rotation of a periodic grid, and it composes cleanly with the multi-level extension in 9.1 (each channel gets its own full-range rank matrix, from which that channel's own cut points are derived). For especially demanding registration/moiré requirements, an explicit cross-channel decorrelation *constraint* can be added to the void-and-cluster construction itself (optimizing not just each channel's own energy map but a joint energy term that penalizes overlap between channels' ON sets at matched densities) — this is a legitimate extension but adds real implementation complexity; **start with independent per-channel seeding and only add joint decorrelation if validation (Section 11.6) shows it is insufficient for the actual target ink set and substrate.** --- ## 10. Calibration and Tone-Response-Curve Integration ### 10.1 Keep the mask device-independent The generated rank matrix / threshold array should be treated as a **device-independent artifact**, versioned and stored separately from any device-, ink-, or substrate-specific calibration data. The recommended pipeline separates two composable stages: ``` source tone value device-specific TRC / dot-gain I(x, y) ──────────────► I'(x,y) = TRC(I(x,y)) ──────────────► compare against T(x,y) ──► ON/OFF (1D lookup, per channel, (blue-noise threshold per device/substrate/ink) array, static, reusable) ``` 1. `T(x, y)` — the blue-noise threshold array from Section 7, generated once, reused across jobs and (where applicable) across devices. 2. A 1D **tone response curve (TRC)** / dot-gain compensation LUT applied to the *input* tone value (equivalently, an inverse curve could instead be applied to the threshold array, but applying it to the input is the more conventional and more easily-audited placement) before the comparison in `apply_screen`. This curve is calibrated per device/substrate/ink combination, typically by printing a ramp, measuring actual optical density or reflectance at each nominal level, and inverting the measured response. This separation is what allows the (computationally nontrivial, `O(N⁴)`) mask-generation step to be a one-time cost, decoupled from the (comparatively cheap, routine) recalibration that happens whenever paper stock, ink batch, or press condition changes. **Do not bake device-specific dot gain into the threshold array itself** — that would require regenerating the mask on every recalibration, throwing away the whole point of precomputing it. ### 10.2 Where TRC calibration itself comes from Full TRC/ICC calibration methodology is outside this document's scope (it belongs to color management, not screening theory), but the interface contract this white paper's deliverable needs to expose is simple: `apply_screen` must accept the *already TRC-corrected* tone value as input, and must not itself apply any nonlinearity beyond the pure threshold comparison. Whatever calibration subsystem exists (or will exist) in the broader RIP pipeline should own the TRC step entirely upstream of this module. --- ## 11. Validation and Testing Methodology A generated mask must pass all of the following before being accepted into production. These are listed in increasing order of subtlety; the first three are cheap correctness checks that should run as automated unit tests on every generation, the later ones are the actual perceptual-quality gate and should run at minimum once per new parameter configuration (new `σ`, new `N`, new `seed_fraction`). ### 11.1 Structural correctness (automated, every run) - **Valid permutation**: every integer `0..N²−1` appears in `rank_matrix` exactly once. (Already asserted in the reference pseudocode, Section 7.3 — keep this assertion in production code, not just during development.) - **Threshold array range**: `T(x,y) ∈ [0, 2^bit_depth − 1]`, and the mapping from rank to threshold is monotonic (guaranteed by construction, but worth a direct test). ### 11.2 Histogram uniformity per threshold slice (automated, every run) For a representative set of gray levels `g ∈ {0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95}`, thresholding a uniform image of value `g` against the tile must produce an ON-pixel count matching `round(g · N²)` exactly (this is guaranteed by construction if the rank-to-threshold mapping is implemented correctly — assert it as a regression test, since a subtle off-by-one in the rank/threshold formula is a realistic implementation bug and this test catches it cheaply). ### 11.3 Radially-averaged power spectrum (RAPS) — the core perceptual check For each representative gray level from 11.2: 1. Threshold a **tiled repetition** of the mask (e.g. 4×4 or 8×8 repetitions of the base tile, to get adequate spectral resolution — a single small tile's own FFT has too few frequency bins to characterize the spectrum meaningfully) at that gray level. 2. Compute the 2D FFT, take the squared magnitude (power spectrum), and radially average it into 1D bins by spatial-frequency radius. 3. Confirm: - **No significant energy below `f_g`** (Section 5.2's formula) — plot or numerically integrate the low-frequency band and confirm it is close to the noise floor. - **A reasonably sharp transition** at `f_g`, not a slow ramp. - **A roughly flat plateau above `f_g`** up to Nyquist — no strong narrow peaks (which would indicate residual periodicity, most likely from too-small a tile relative to the repetition count used in step 1, or from Bayer-like artifacts in the underlying rank matrix if the algorithm was implemented incorrectly). This is the standard blue-noise sanity check used throughout the source literature (Ulichney 1988, 1993) and is the single most diagnostic test in this checklist — most implementation defects (missing toroidal wrap, wrong phase logic, wrong σ) show up clearly here even when the mask "looks fine" on casual visual inspection. ### 11.4 Anisotropy check (automated, at least once per parameter configuration) Compute the **full 2D** power spectrum (not just its radial average) for the mid-tone (`g = 0.5`) case, and compare energy in angular sectors (e.g. eight 45°-wide sectors) at matched radius. A well-formed blue-noise mask should show comparable energy across sectors; a systematic excess along the horizontal/vertical or diagonal directions indicates directional bias (Section 5.4) — commonly traceable to a non-circular convolution kernel, an incorrectly-implemented toroidal wrap, or (less commonly) an anisotropic RNG artifact in Phase 0's initial pattern. ### 11.5 Tileability / seam check (visual, at least once per parameter configuration) Render a 3×3 (or larger) grid of the tile at a mid gray level and visually and numerically (e.g. comparing pixel statistics in a band straddling each internal tile boundary against statistics from equivalent interior regions) confirm no seam, edge artifact, or periodic discontinuity is visible or measurable at tile boundaries. ### 11.6 Ramp test (visual, at least once per parameter configuration) Render a full 0→100% gradient ramp at production resolution through `apply_screen` and inspect for banding, false contouring, or abrupt texture/character changes between adjacent gray levels — this is the most direct end-to-end check that the three-phase construction (Section 6.4) actually delivered a smooth, monotonic quality progression across the whole tonal range, not just at the levels explicitly checked in 11.2–11.3. ### 11.7 Cross-channel test (multi-channel deployments only) Per Section 9.2: overlay two or more independently-seeded channel masks at matched densities (e.g. both at 50%) and check for visible moiré, coincident dot alignment, or elevated combined graininess relative to either channel alone. Automate this as a combined RAPS/anisotropy check on the union pattern where practical, backed by a visual spot-check. --- ## 12. Production Integration Notes - **Precompute once, ship as a static asset.** Generate the rank matrix (and/or the derived threshold array) offline; do not regenerate at runtime. Store the canonical rank matrix (not just a derived threshold array at one bit depth) so higher- or lower-precision threshold arrays can be re-derived later without rerunning the `O(N⁴)` generation step. - **Storage format.** A raw binary array (row-major, fixed-width integers) or a lossless image format (16-bit PNG, if 16-bit precision is needed) are both reasonable; avoid any lossy format for the threshold array itself. - **Runtime cost.** `apply_screen`, per Section 7.3, is a single modulo-indexed threshold comparison per pixel — trivially parallelizable and GPU-friendly (a natural fit for a Vulkan/compute-shader RIP pipeline: the threshold tile is a small, cacheable texture sampled with wrap addressing, and the comparison itself is a single scalar operation per fragment/thread). - **Version and log generation parameters.** Store `(size, sigma, seed_fraction, rng_seed, algorithm_version)` alongside the shipped mask, both for reproducibility and so that a future regression (e.g. after a dependency update changes floating-point behavior subtly) can be diagnosed by regenerating from the same parameters and diffing against the archived rank matrix. - **Do not couple mask generation to per-job calibration.** Per Section 10, keep TRC/dot-gain entirely upstream of `apply_screen`; a mask regeneration should only be triggered by an actual change to `(size, sigma, seed_fraction)` — a deliberate quality/design decision — never by routine recalibration. --- ## 13. Glossary - **Area coverage** — the fraction of a local neighborhood covered by marks; the physical quantity that halftoning uses to approximate a continuous tone value. - **Dither array / threshold array** — a precomputed matrix of threshold values, tiled across an image, used for ordered dithering. - **Dither pattern / rank matrix** — the underlying combinatorial ordering (0..N²−1) from which a threshold array is derived by rescaling. - **RAPS** — radially-averaged power spectrum; the standard 1D characterization of a 2D pattern's frequency content, used as the primary blue-noise diagnostic. - **Principal frequency (`f_g`)** — the density-dependent spatial frequency at which a blue-noise pattern's spectral energy is concentrated; see Section 5.2. - **Toroidal / circular convolution** — convolution with wrap-around boundary conditions, required for seamless tileability. - **Void-and-cluster** — the spatial-domain, greedy algorithm (Ulichney 1993) specified in Section 6 for constructing a blue-noise rank matrix. - **Direct Binary Search (DBS)** — a higher-cost, iterative-optimization alternative to void-and-cluster; out of scope for this specification (Section 6.6). - **Tone response curve (TRC)** — the device/substrate-specific nonlinear mapping between nominal and actual (measured) tone, applied upstream of screening (Section 10). --- ## 14. References - Floyd, R. W., Steinberg, L. (1976). *An Adaptive Algorithm for Spatial Grayscale.* Proc. SID, 17, 75–77. - Jarvis, J. F., Judice, C. N., Ninke, W. H. (1976). *A Survey of Techniques for the Display of Continuous Tone Pictures on Bilevel Displays.* Computer Graphics and Image Processing, 5, 13–40. - Knuth, D. E. (1987). *Digital Halftones by Dot Diffusion.* ACM Transactions on Graphics, 6(4), 245–273. - Ulichney, R. A. (1987). *Digital Halftoning.* MIT Press. - Ulichney, R. A. (1988). *Dithering with Blue Noise.* Proceedings of the IEEE, 76(1), 56–79. - Ulichney, R. A. (1993). *The Void-and-Cluster Method for Dither Array Generation.* Proc. SPIE 1913, Human Vision, Visual Processing, and Digital Display IV, 332–343. - Ulichney, R. A. (2000). *A Review of Halftoning Techniques.* Proc. SPIE 3963. - Mitsa, T., Parker, K. J. (1992). *Digital Halftoning Technique Using a Blue-Noise Mask.* Journal of the Optical Society of America A, 9(11). - Ostromoukhov, V. (2001). *A Simple and Efficient Error-Diffusion Algorithm.* Proc. SIGGRAPH 2001. - Ostromoukhov, V. (1996). *Multi-Level Colour Halftoning Algorithms.* LIRIS / CNRS. - Kang, H. R. (1999). *Digital Color Halftoning.* SPIE Press / Wiley. - Allebach, J. P., and collaborators — various papers on Direct Binary Search halftoning (model-based, clustered-dot, and multiscale variants), cited comprehensively in the modern literature review below. - "Digital Halftoning: A Literature Review from Classical Methods to Deep Reinforcement Learning" (2024/2025) — comprehensive modern survey covering DBS variants and reinforcement-learning-based screening; recommended further reading beyond this document's scope. --- ## Appendix A: Reference Pseudocode Summary The full generation pseudocode is given inline in Section 7.3; this appendix summarizes the call graph for quick reference during implementation review: ``` void_and_cluster_rank(size, sigma, seed_fraction, rng_seed) ├─ make_gaussian_patch(sigma) ├─ Phase 0: build prototype │ ├─ random initial pattern (n0 ON pixels) │ ├─ toggle_pixel(...) × n0 # build initial energy map │ └─ loop: tightest_cluster / largest_void / toggle_pixel # until converged ├─ Phase 1: rank 0 .. n0-1 │ └─ loop: tightest_cluster / toggle_pixel (OFF) # counting down ├─ Phase 2: rank n0 .. N²/2-1 │ └─ loop: largest_void / toggle_pixel (ON) # counting up ├─ Phase 3: rank N²/2 .. N²-1 │ ├─ build complement pattern + fresh energy map │ └─ loop: tightest_cluster (on complement) / toggle_pixel (OFF, on complement) └─ assert valid permutation rank_to_threshold(rank_matrix, bit_depth) apply_screen(image, threshold_tile) validate_blue_noise(threshold_tile, levels) ├─ structural correctness checks (11.1) ├─ histogram uniformity (11.2) ├─ RAPS computation + low-frequency / flatness checks (11.3) ├─ anisotropy sector comparison (11.4) └─ (tileability / ramp / cross-channel checks are integration-level, not unit-level) ``` ## Appendix B: Suggested Test Suite Structure ``` tests/ test_prototype_convergence.py # Phase 0 terminates within max_iters, for a range of (N, sigma, seed_fraction) test_rank_matrix_permutation.py # 11.1 test_histogram_uniformity.py # 11.2, across the representative gray-level set test_raps_blue_noise.py # 11.3, asserts low-freq energy below a tolerance, plateau flatness within a tolerance test_anisotropy.py # 11.4, sector-energy comparison within a tolerance test_tileability_seam.py # 11.5, statistical comparison across tile boundaries test_ramp_no_banding.py # 11.6, automated proxy (e.g. local variance continuity) plus a saved visual artifact for manual review test_multichannel_decorrelation.py # 11.7, for multi-channel deployments test_determinism.py # 7.5, same params + seed -> byte-identical rank matrix test_apply_screen_tiling.py # apply_screen correctness on non-tile-aligned image sizes ``` *End of white paper.* --- © Hybrid Software Helix. Confidential document — unauthorized reproduction or distribution prohibited.