{'fraction_non_alphanumeric': 0.03534717494894486, 'fraction_numerical': 0.007981620149761743, 'mean_word_length': 4.6655739285542115, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perceptual error optimization for Monte Carlo rendering

Synthesizing realistic images involves computing high-dimensional light-transport integrals. In practice, these integrals are numerically estimated via Monte Carlo integration. The error of this estimation manifests itself as conspicuous aliasing or noise. To ameliorate such artifacts and improve image fidelity, we propose a perception-oriented framework to optimize the error of Monte Carlo rendering. We leverage models based on human perception from the halftoning literature. The result is an optimization problem whose solution distributes the error as visually pleasing blue noise in image space. To find solutions, we present a set of algorithms that provide varying trade-offs between quality and speed, showing substantial improvements over prior state of the art. We perform evaluations using quantitative and error metrics, and provide extensive supplemental material to demonstrate the perceptual improvements achieved by our methods.

. We devise a perceptually based model to optimize the error of Monte Carlo renderings. Here we show our vertical iterative minimization algorithm from Section 4.1: Given 4 input samples per pixel (spp), it selects a subset of them to produce an image with substantially improved visual fidelity over a simple 4-spp average. The optimization is guided by a surrogate image obtained by regularizing the noisy input; we also show using the ground-truth image as a guide. The power spectrum of the image error, computed on 32×32-pixel tiles, indicates that our method distributes pixel error with locally blue-noise characteristics.
Synthesizing realistic images involves computing high-dimensional lighttransport integrals. In practice, these integrals are numerically estimated via Monte Carlo integration. The error of this estimation manifests itself as conspicuous aliasing or noise. To ameliorate such artifacts and improve image fidelity, we propose a perception-oriented framework to optimize the error of Monte Carlo rendering. We leverage models based on human perception from the halftoning literature. The result is an optimization problem whose solution distributes the error as visually pleasing blue noise in image space. To find solutions, we present a set of algorithms that provide varying trade-offs between quality and speed, showing substantial improvements over prior state of the art. We perform evaluations using quantitative and error metrics, and provide extensive supplemental material to demonstrate the perceptual improvements achieved by our methods.

INTRODUCTION
Monte Carlo sampling produces approximation error. In rendering, this error can cause visually displeasing image artifacts, unless control is exerted over the correlation of the individual pixel estimates. A standard approach is to decorrelate these estimates by randomizing the samples independently for every pixel, turning potential structured artifacts into white noise.
In digital halftoning, the error induced by quantizing continuoustone images has been studied extensively. Such studies have shown that a blue-noise distribution of the quantization error is perceptually optimal [Ulichney 1987], achieving substantially higher image fidelity than a white-noise distribution. Recent works have proposed empirical means to transfer these ideas to image synthesis Ahmed and Wonka 2020]. Instead of randomizing the pixel estimates, these methods introduce negative correlation between neighboring pixels, exploiting the local smoothness in images to push the estimation error to the high-frequency spectral range.
We propose a theoretical formulation of perceptual error for image synthesis which unifies prior methods in a common framework and formally justifies the desire for blue-noise error distribution. We extend the comparatively simpler problem of digital halftoning [Lau and Arce 2007] where the ground-truth image is given, to the substantially more complex one of rendering where the ground truth is the sought result and thus unavailable. Our formulation bridges the gap between multi-tone halftoning and rendering by interpreting Monte Carlo estimates for a pixel as its admissible 'quantization levels'. This insight allows virtually any halftoning method to be adapted to rendering. We demonstrate this for the three main classes of halftoning algorithms: dither-mask halftoning, error diffusion halftoning, and iterative energy minimization halftoning.
Existing methods  can be seen as variants of dither-mask halftoning. They distribute pixel error according to masks that are optimized w.r.t. a target kernel, typically a Gaussian. The kernel can be interpreted as an approximation to the human visual system's point spread function [Daly 1987;Pappas and Neuhoff 1999]. We revisit the kernel-based perceptual model from halftoning [Sullivan et al. 1991;Pappas and Neuhoff 1999] and adapt it to rendering. The resulting energy can be directly used for optimizing Monte Carlo error distribution without the need for a mask. This formulation help us expose the underlying assumptions of existing methods and quantify their limitations. In summary: • We formulate an optimization problem for rendering error by leveraging kernel-based perceptual models from halftoning. • Our formulation unifies prior blue-noise error distribution methods and makes all their assumptions explicit, outlining general guidelines for devising new methods in a principled manner.
• Unlike prior methods, our formulation simultaneously optimizes for both the magnitude and the image distribution of pixel error.
• We devise four different practical algorithms based on iterative minimization, error diffusion, and dithering from halftoning.
• We demonstrate substantial visual improvements over prior art, while using the same input rendering data.

RELATED WORK
Our work focuses on reducing and optimizing the distribution of Monte Carlo pixel-estimation error. In this section we review prior work with similar goals in digital halftoning (Section 2.1) and image synthesis guided by energy-based (Section 2.2) and perception-based (Section 2.3) error metrics. We achieve error reduction through careful sample placement and processing, and discuss related rendering approaches (Section 2.4).

Digital halftoning
Digital halftoning [Lau and Arce 2007] involves creating the illusion of continuous-tone images through the arrangement of binary elements; various algorithms target different display devices. Bayer [1973] developed the widely used dispersed-dot ordered-dither patterns. Allebach and Liu [1976] introduced the use of randomness in clustered-dot ordered dithering. Ulichney [1987] introduced bluenoise patterns that yield better perceptual quality, and Mitsa and Parker [1991] mimicked those patterns to produce dither arrays (i.e., masks) with high-frequency characteristics. Sullivan et al. [1991] developed a Fourier-domain energy function to obtain visually optimal halftone patterns; the optimality is defined w.r.t. computational models of the human visual system.  devised a practical algorithm for blue-noise dithering through a spatial-domain interpretation of Sullivan et al.'s model. Their approach was later refined by Pappas and Neuhoff [1999].
The void-and-cluster algorithm [Ulichney 1993] uses a Gaussian kernel to create dither masks with isotropic blue-noise distribution. This approach has motivated various structure-aware halftoning algorithms in graphics [Ostromoukhov 2001;Pang et al. 2008;Chang et al. 2009]. In the present work, we leverage the kernel-based model Pappas and Neuhoff 1999] in the context of Monte Carlo rendering ].

Quantitative error assessment in rendering
It is convenient to measure the error of a rendered image as a single value; vector norms like the mean squared error (MSE) are most commonly used. However, it is widely acknowledged that such simple metrics do not accurately reflect visual quality as they ignore the perceptually important spatial arrangement of pixels. Various theoretical frameworks have been developed in the spatial [Niederreiter 1992;Kuipers and Niederreiter 1974] and Fourier [Singh et al. 2019] domains to understand the error reported through these metrics. The error spectrum ensemble [Celarek et al. 2019] measures the frequency-space distribution of the error.
Many denoising methods [Zwicker et al. 2015] employ the aforementioned metrics to obtain noise-free results from noisy renderings. Even if the most advanced denoising techniques driven by such metrics can efficiently steer adaptive sampling [Chaitanya et al. 2017;Kuznetsov et al. 2018;Kaplanyan et al. 2019], they locally determine the number of samples per pixel, ignoring the aspect of their specific layout in screen space.
Our optimization framework employs a perceptual MSE-based metric that accounts for both the magnitude and the spatial distribution of pixel-estimation error. We argue that the spatial sample layout plays a crucial role in the perception of a rendered image; the most commonly used error metrics do not capture this aspect.

Perceptual error assessment in rendering
The study of the human visual system (HVS) is still ongoing, and well understood are mostly the early stages of the visual pathways from the eye optics, through the retina, to the visual cortex. This limits the scope of existing HVS computational models used in imaging and graphics. Such models should additionally be computationally efficient and generalize over the simplistic stimuli that have been used in their derivation through psychophysical experiments.
Contrast sensitivity function. The contrast sensitivity function (CSF) is one of the core HVS models that fulfills the above conditions and comprehensively characterizes overall optical [Westheimer 1986;Deeley et al. 1991] and neural [Souza et al. 2011] processes in detecting contrast visibility as a function of spatial frequency. While originally modeled as a band-pass filter [Barten 1999;Daly 1992], the CSF's shape changes towards a low-pass filter with retinal eccentricity [Robson and Graham 1981;Peli et al. 1991] and reduced luminance adaptation in scotopic and mesopic levels [Wuerger et al. 2020]. Low-pass characteristics are also inherent for chromatic CSFs [Mullen 1985;Wuerger et al. 2020;Bolin and Meyer 1998]. In many practical imaging applications, e.g., JPEG compression [Rashid et al. 2005], rendering [Ramasubramanian et al. 1999], or halftoning [Pappas andNeuhoff 1999], the CSF is modeled as a low-pass filter, which also allows for better control of image intensity. By normalizing such a CSF by the maximum contrastsensitivity value, a unitless function akin to the modulation transfer function (MTF) can be derived [Daly 1987;Mannos and Sakrison 1974;Mantiuk et al. 2005;Sullivan et al. 1991;Souza et al. 2011] that after transforming from the frequency to the spatial domain results in the point spread function (PSF) Pappas and Neuhoff 1999]. Following Pappas and Neuhoff [1999], we approximate such a PSF by a Gaussian filter; the resulting error is practically negligible for a pixel density of 300 dots per inch (dpi) and observer-to-screen distance larger than 60 cm.
Advanced quality metrics. More costly, and often less robust, modeling of the HVS beyond the CSF is performed in advanced quality metrics [Lubin 1995;Daly 1992;Mantiuk et al. 2011]. Such metrics have been adapted to rendering to guide the computation to image regions where the visual error is most strongly perceived Meyer 1995, 1998;Ramasubramanian et al. 1999;Ferwerda et al. 1996;Myszkowski 1998;Volevich et al. 2000]. An important application is visible noise reduction in path tracing via content-adaptive sample-density control Meyer 1995, 1998;Ramasubramanian et al. 1999]. Our framework enables significant reduction of noise visibility for the same sampling budget. Mitchell [1991] first observed that high-frequency error distribution is desirable for stochastic rendering. Only recently,  adopted techniques from halftoning to correlate pixel samples in screen space and distribute path-tracing error as blue noise, with substantial perceptual quality improvements.  built on this idea to develop a progressive quasi-Monte Carlo sampler that further improves quality. Ahmed and Wonka [2020] proposed a technique to coordinate quasi-Monte Carlo samples in screen space inspired by error diffusion.

Blue-noise error distribution in rendering
Motivated by the results of ,  devised a method to directly optimize the distribution of pixel estimates, without operating on individual samples. Their pixel permutation strategy fits the initially white-noise pixel intensities to a prescribed blue-noise mask. This approach scales well with sample count and dimension, though its reliance on prior pixel estimates makes it practical only for animation rendering where it is susceptible to quality degradation.
We propose a perceptual error framework that unifies these two general approaches, exposing the assumptions of existing methods and providing guidelines to alleviate some of their drawbacks.

PERCEPTUAL ERROR MODEL
Our aim is to produce Monte Carlo renderings that, at a fixed sampling rate, are perceptually as close to the ground truth as possible. This goal requires formalizing the perceptual image error along with an optimization problem that minimizes it. In this section, we build a perceptual model upon the extensive studies done in the halftoning literature. We will discuss how to efficiently solve the resulting optimization problem in Section 4.
Given a ground-truth image and its quantized or stochastic approximation , we denote the (signed) error image by To minimize the error, it is convenient to quantify it as a single value. A common approach is to take the ℒ 1 , ℒ 2 , or ℒ∞ norm of the Image Image spectrum Kernel spectrum Product spectrum Error images w and b with respective white-noise, ⋃︀ˆw⋃︀ 2 , and bluenoise, ⋃︀ˆb⋃︀ 2 , power spectra. For a low-pass kernel modeling the PSF of the HVS (here a Gaussian with std. dev. = 1), the product of its spectrum ⋃︀ˆ⋃︀ 2 with ⋃︀ˆb⋃︀ 2 has lower magnitude than the product with ⋃︀ˆw⋃︀ 2 . This corresponds to lower perceptual sensitivity to b , even though w has the same amplitude as it is obtained by randomly permuting the pixels of b .
image interpreted as a vector. Such simple metrics are permutationinvariant, i.e., they account for the magnitudes of individual pixel errors but not for their distribution over the image. This distribution is an important factor for the perceived fidelity, since contrast perception is an inherently spatial characteristic of the HVS (Section 2.3). Our model is based on perceptual halftoning metrics that capture both the magnitude and the distribution of error.

Motivation
Halftoning metrics model the processing done by the HVS as a convolution of the error image with a kernel : The convolution is equivalent to the element-wise product of the corresponding Fourier spectraˆandˆ, whose 2-norm in turn equals the inner product of the power spectra images ⋃︀ˆ⋃︀ 2 and ⋃︀ˆ⋃︀ 2 . Sullivan et al. [1991] optimized the error image to minimize the error (2) w.r.t. a kernel that approximates the HVS's modulation transfer function ⋃︀ˆ⋃︀ (MTF) [Daly 1987].  used a similar model in the spatial domain with a kernel that approximates the PSF 1 of the human eye. That kernel is low-pass, and the optimization naturally yields blue-noise 2 distribution in the error image , as we show later in Fig. 5. The blue-noise distribution can thus be seen as byproduct of the optimization which pushes the spectral components of the error to the frequencies least visible to the human eye (see Fig. 2). To better understand the spatial aspects of contrast sensitivity in the HVS, the MTF is usually modeled over a range of viewing distances [Daly 1992]. This is done to account for the fact that with increasing viewer distance, spatial frequencies in the image are = 0 = 0.25 = 0.5 = 1 Fig. 3. The appearance of blue noise (left images) converges to a constant image faster than white noise (right images) with increasing observer distance, here emulated via the standard deviation of a Gaussian kernel. We provide a formal connection between and the viewing distance in Section 6.
projected to higher spatial frequencies onto the retina. These frequencies eventually become invisible, filtered out by the PSF which expands its corresponding kernel in image space. We recreate this experiment to see the impact of distance on the image error. In Fig. 3, we convolve white-and blue-noise distributions with a Gaussian kernel of increasing standard deviation corresponding to increasing observer-to-screen distance. The high-frequency blue-noise distribution reaches a homogeneous state (where the tone appears constant) faster compared to the all-frequency white noise. This means that high-frequency error becomes indiscernible at closer viewing distances, where the HVS ideally has not yet started filtering out actual image detail which is typically low-to mid-frequency. In Section 6 we discuss how the kernel's standard deviation encodes the viewing distance w.r.t. to the screen resolution.

Our model
In rendering, the value of each pixel is a light-transport integral. Point-sampling its integrand with a sample set yields a pixel estimate ( ). The signed pixel error is thus a function of the sample set: ( ) = ( ) − , where is the reference (i.e., groundtruth) pixel value. The error of the entire image can be written as where = { 1 , . . . , } is an "image" containing the sample set for all pixels. With these definitions, we can express the perceptual error in Eq. (2) for the case of Monte Carlo rendering as a function of the sample-set image , given a kernel : Our goal is to minimize the perceptual error (4). We formulate this task as an optimization problem: The minimizing sample-set image yields an image estimate ( ) that is closest to the reference w.r.t. the kernel . The search space Ω Ω Ω is the set of all possible locations for every sample of every pixel. The total number of samples in is typically bounded by a given target sampling budget. Practical considerations may also restrict the search space Ω Ω Ω, as we will exemplify in the following section. Note that the classical MSE metric corresponds to using a zerowidth (i.e., one-pixel) kernel in Eq. (4). However, the MSE accounts only for the magnitude of the error , while using wider kernels (such as the PSF) accounts for both magnitude and distribution. Consequently, while the MSE can be minimized by optimizing pixels independently, minimizing the perceptual error requires coordination between pixels. In the following section, we devise strategies for solving this optimization problem.

DISCRETE OPTIMIZATION
In our optimization problem (5), the search space for each sample in every pixel is a high-dimensional unit hypercube. Every point in this so-called primary sample space maps to a light-transport path in the scene [Pharr et al. 2016]. Optimizing for the sample-set image thus entails evaluating the contributions ( ) of all corresponding paths. This evaluation is costly, and for any non-trivial scene, is a function with complex shape and many discontinuities. This precludes us from studying all (uncountably infinite) sample locations in practice.
To make the problem tractable, we restrict the search in each pixel to a finite number of (pre-defined) sample sets. We devise two variants of the resulting discrete optimization problem, which differ in their definition of the global search space Ω Ω Ω. In the first variant, each pixel has a separate list of sample sets to choose from ("vertical" search space). The setting is similar to that of (multitone) halftoning [Lau and Arce 2007], which allows us to import classical optimization techniques from that field, such as iterative minimization, error diffusion, and mask-based dithering. In the second variant, each pixel has one associated sample set, and the search space comprises permutations of these assignments ("horizontal" search space). We develop a greedy iterative optimization method for this second variant.
In contrast to halftoning, in our setting the ground-truth image -required to compute the error image during optimization-is not readily available. Below we describe our algorithms assuming the ground truth is available; in Section 5 we will discuss how to substitute it with a surrogate to make the algorithms practical.

Vertical search space
Our first variant considers a "vertical" search space where the sample set for each of the image pixels is one of given sets: 3 The objective is to find a sample set for every pixel such that all resulting pixel estimates together minimize the perceptual error (4).
This problem is almost identical to that of multi-tone halftoning. The difference is that in our setting the "quantization levels", i.e., the pixel estimates, are distributed non-uniformly and vary per pixel as they are not fixed but are the result of point-sampling a light-transport integral. This similarity allows us to directly apply existing optimization techniques from halftoning. We consider three such methods, which we outline in Alg. 1 and describe next.
Iterative minimization. State-of-the-art halftoning methods attack the problem (7) directly via greedy iterative minimization Pappas and Neuhoff 1999]. After initializing every pixel to a random quantization level, we traverse the image in serpentine order (as is standard practice in halftoning) and for each pixel choose the level that minimizes the energy. Several full-image iterations are performed; in our experiments convergence to a local minimum is achieved within 10-20 iterations.
As a further improvement, the optimization can be terminated when no pixels are updated within one full iteration, or when the perceptual-error reduction rate drops below a certain threshold. Traversing the pixels in random order allows terminating at any point but converges slightly slower.
Error diffusion. A classical halftoning algorithm, error diffusion scans the image pixel by pixel, snapping each reference value to the closest quantization level and distributing the resulting pixel error to yet-unprocessed nearby pixels according to a given kernel . We use the empirically derived kernel of Floyd and Steinberg [1976] which has been shown to produce an output that approximately minimizes Eq. (7) [Hocevar and Niger 2008]. Error diffusion is faster than iterative minimization but yields less optimal solutions.
Dithering. The fastest halftoning approach quantizes pixel values using thresholds stored in a pre-computed dither mask (or matrix) [Spaulding et al. 1997]. For each pixel, the two quantization levels that tightly envelop the reference value (in terms of brightness) are found, and one of the two is chosen based on the threshold assigned to the pixel by the mask.
Dithering can be understood as performing the perceptual error minimization in two steps. First, an offline optimization encodes the error distribution optimal for the target kernel into a mask. Then, for a given image, the error magnitude is minimized by restricting the quantization to the two closest levels per pixel, and the maskdriven choice between them applies the target distribution of error.

Horizontal search space
We now describe the second, "horizontal" discrete variant of our minimization formulation (5). It considers a single sample set assigned to each of the pixels, all represented together as a sampleset image . The search space comprises all possible permutations Π( ) of these assignments: The goal is to find a permutation ( ) that minimizes the perceptual error (4). The optimization problem (5) thus takes the form min ∈ Π( ) ∏︁ * ( ( ( )) − )∏︁ 2 2 .
Algorithm 1. Three algorithms to (approximately) solve the vertical search space optimization problem (7). The output is an image = { 1 , . . . , }, given a reference image and a stack of initial image estimates 1 , . . . , . Iterative minimization updates pixels repeatedly, for each selecting the estimate that minimizes the perceptual error (4). Error diffusion quantizes each pixel to the closest estimate, distributing the error to its neighbors based on a kernel . Dithering quantizes each pixel in based on thresholds looked up in a dither mask (optimized w.r.t. the kernel ).
1: function IterativeMinimization( , , 1 , . . . , We can explore the permutation space Π( ) by swapping the sampleset assignments between pixels. The minimization requires updating the image estimate ( ( )) for each permutation ( ), i.e., after every swap. Such updates are costly as they involve re-sampling both pixels in each of potentially millions of swaps. We need to eliminate these extra rendering invocations during the optimization to make it practical. To that end, we observe that for pixels solving similar light-transport integrals, swapping their sample sets gives a similar result to swapping their estimates. We therefore restrict the search space to permutations that can be generated through swaps between such (similar) pixels. This enables an efficient optimization scheme that directly swaps the pixel estimates of an initial rendering ( ).
Error decomposition. Formally, we express the estimate produced by a sample-set permutation in terms of permuting the pixels of the initial rendering: ( ( )) = ( ( )) + Δ Δ Δ( ). The error Δ Δ Δ is zero when the swapped pixels solve the same integral. Substituting into Eq. (9), we can approximate the perceptual error by (see Appendix A) ( ) = ∏︁ * ( ( ( )) − + Δ Δ Δ( ))∏︁ 2 2 (10a) ≈ ∏︁ * ( ( ( )) − )∏︁ 2 2 + ∏︁ ∏︁ 2 1 ∑ ( , ( )) = ( ), (10b) where we write the error ( ) as a function of only, to emphasize that everything else is fixed during the optimization. In the approximation , the term ( , ( )) measures the dissimilarity between pixel and the pixel ( ) it is relocated to by the permutation. The purpose of this metric is to predict how different we expect the result of re-estimating the pixels after swapping their sample sets to be compared to simply swapping their initial estimates. It can be constructed based on knowledge or assumptions about the image.
Local similarity assumption. Our implementation uses a simple binary dissimilarity function that returns zero when and ( ) are within some distance and infinity otherwise. We set ∈ (︀1, 3⌋︀; it should ideally be locally adapted to the image smoothness. This allows us to restrict the search space Π( ) only to permutations that swap adjacent pixels where it is more likely that Δ Δ Δ is small. More elaborate heuristics could better account for pixel (dis)similarity.
Iterative minimization. We devise a greedy iterative minimization scheme for this horizontal formulation, similar to the one in Alg. 1. Given an initial image estimate ( ), produced by randomly assigning a sample set to every pixel, our algorithm goes over all pixels and for each considers swaps within a (2 + 1) 2 neighborhood; we use = 1. The swap that brings the largest reduction in the perceptual error is accepted. Algorithm 2 provides pseudocode. In our experiments we run = 10 full-image iterations. As before, the algorithm could be terminated based on the swap reduction rate or the error reduction rate. We explore additional optimizations in supplemental Section 3.
The parameter balances between the cost of one iteration and the amount of exploration it can do. Note that this parameter is different from the maximal relocation distance in the dissimilarity metric, with ≤ .
Due to the pixel (dis)similarity assumptions, the optimization can produce some mispredictions, i.e., it may swap the estimates of pixels for which swapping the sample sets produces a significantly different result. Thus the image ( ( )) cannot be used directly as a final estimate. We therefore re-render the image using the optimized permutation to obtain the final estimate ( ( )).

Discussion
Search space. We discretize the search space Ω Ω Ω to make the optimization problem (5) tractable. To make it truly practical, it is also necessary to avoid repeated image estimation (i.e., ( ) evaluation) during the search for the solution . Our vertical (7) and horizontal (9) optimization variants are formulated specifically with this goal in mind. All methods in Algs. 1 and 2 operate on pre-generated image estimates that constitute the solution search space.
Our vertical formulation takes a collection of input estimates { , = ( , )} =1 for every pixel , one for each sample set , .
Noting that , are MC estimates of the true pixel value, this collection can be cheaply expanded to a size as large as 2 − 1 by Algorithm 2. Given a convolution kernel , a reference image , an initial sample-set assignment , and an image estimate ( ) computed with that assignment, our greedy algorithm iteratively swaps sample-set assignments between neighboring pixels to minimize the perceptual error (10b), producing a permutation of the initial assignment. if ( ⇆ ( )) < ( ′ ( )) then = ′ taking the average of the estimates in each of its subsets (excluding the empty subset). In practice only a fraction of these subsets can be used, since the size of the power set grows exponentially with . It may seem that this approach ends up wastefully throwing away most input estimates. But note that these estimates actively participate in the optimization and provide the space of possible solutions. Carefully selecting a subset per pixel can yield a higherfidelity result than blindly averaging all available estimates, as we will show repeatedly in Section 7.
In contrast, our horizontal formulation builds a search space given just a single input estimate per pixel. We consciously restrict the space to permutations between nearby pixels, so as to leverage local pixel similarity and avoid repeated pixel evaluation during optimization. The disadvantage of this approach is that it requires rerendering the image after optimization, with uncertain results (due to mispredictions) that can lead to local degradation of image quality. Mispredictions can be reduced by exploiting knowledge about the rendering function ( ), e.g., through depth, normal, or texture buffers; we explore this in supplemental Section 2. Additionally, while methods like iterative minimization (Alg. 2) and dithering (Section 5.2) can be adapted to this search space, reformulating other halftoning algorithms such as error diffusion is non-trivial.
A hybrid formulation is also conceivable, taking a single input estimate per pixel (like horizontal methods) and considering a separate (vertical) search space for each pixel constructed by borrowing estimates from neighboring pixels. Such an approach could benefit from advanced halftoning optimization methods, but could also suffer from mispredictions and require re-rendering. We leave the exploration of this approach to future work.
Finally, it is worth noting that discretization is not the only route to practicality. Equation (5) can be optimized on the continuous space Ω Ω Ω if some cheap-to-evaluate proxy for the rendering function is available. Such a continuous approximation may be analytical (based on prior knowledge or assumptions) or obtained by reconstructing a point-wise evaluation. However, continuous-space optimization can be difficult in high dimensions (e.g., number of light bounces) where non-linearities and non-convexity are exacerbated.
Optimization strategy. Another important choice is the optimization method. For the vertical formulation, iterative minimization provides the best flexibility and quality but is the most computationally expensive. Error diffusion and dithering are faster but only approximately solve Eq. (7).
One difference between classical halftoning and our vertical setting is that quantization levels are non-uniformly distributed and differ between pixels. This further increases the gap in quality between the image-adaptive iterative minimization and error diffusion (which can correct for these differences) and the non-adaptive dithering, compared to the halftoning setting. The main advantage of dithering is that it involves the kernel explicitly, while the error-diffusion kernel cannot be related directly to .

PRACTICAL APPLICATION
We now turn to the practical use of our error optimization framework. In both our discrete formulations from Section 4, the search space is determined by a given collection of sample sets , for every pixel , with = 1... (in the horizontal setting = 1). The optimization is then driven by the corresponding estimates , . We consider two ways to obtain these estimates, leading to different practical trade-offs: (1) direct evaluation of the samples by rendering a given scene and (2) using a proxy for the rendering function. We show how prior works correspond to using either approach within our framework, which helps expose their implicit assumptions.

Surrogate for ground truth
The goal of our optimization is to perceptually match an image estimate to the ground truth as closely as possible. Unfortunately, the ground truth is unknown in our setting, unlike in halftoning. The best we can do is substitute it with a surrogate image ′ ′ ′ . Such an image can be obtained either from available pixel estimates or by making assumptions about the ground truth. We will discuss specific approaches in the following, but it is already worth noting that all existing error-distribution methods rely on such a surrogate, whether explicitly or implicitly. And since the surrogate guides the optimization, its fidelity directly impacts the fidelity of the output.

A-posteriori optimization
Given a scene and a viewpoint, initial pixel estimates can be obtained by invoking the renderer with the input samples: , = ( , ). A surrogate can then be constructed from those estimates; in our implementation we denoise the estimate-average image (Section 7.1). Having the estimates and the surrogate, we can run any of the methods in Algs. 1 and 2. Vertical algorithms directly output an image ; horizontal optimization yields a sample-set image that requires an additional rendering invocation: = ( ).
This general approach of utilizing sampled image information was coined a-posteriori optimization by ; they proposed two such methods. Their histogram method operates in a vertical setting, choosing one of the (sorted) estimates for each pixel based on the respective value in a given blue-noise dither mask. Such sampling corresponds to using an implicit surrogate that is the median estimate for every pixel, which is what the mean of the dither mask maps to. Importantly, any one of the estimates for a pixel can be selected, whereas in classical dithering the choice is between the two quantization levels that tightly envelop the reference value (Section 4.1) [Spaulding et al. 1997]. Such selection can yield large error, especially for pixels whose corresponding mask values deviate strongly from the mask mean. This produces image fireflies that do not appear if simple estimate averages are taken instead (see Fig. 6).
The permutation method of  operates in a horizontal setting. Given an image estimate, it finds pixel permutations within small tiles that minimize the distance between the estimates and the values of a target blue-noise mask. This matching transfers the mask's distribution to the image signal rather than to its error. The two are equivalent only when the signal within each tile is constant. The implicit surrogate in this method is thus a tile-wise constant image (shown more formally in supplemental Section 5). In our framework the use of a surrogate is explicit, which enables full control over the quality of the error distribution.

A-priori optimization
Optimizing perceptual error is possible even in the absence of information about a specific image. In our framework, the goal of such an a-priori approach (as coined by ) is to compute a sample-set image by using surrogates for both the ground-truth image and the rendering function ( ), constructed based on smoothness assumptions. The samples can then produce a high-fidelity estimate of any image that meets those assumptions.
Lacking prior knowledge, one could postulate that every pixel has the same rendering function: (⋅) = (⋅); the image surrogate ′ ′ ′ is thus constant. While in practice this assumption (approximately) holds only locally, the optimization kernel is also expected to have compact support. The shape of can be assumed to be (piecewise) smooth and approximable by a cheap analytical function ′ .
With the above surrogates in place, we can run our algorithms to optimize a sample-set image . The constant-image assumption makes horizontal algorithms well-suited for this setting as it makes the swapping-error term Δ Δ Δ in Eq. (10a) vanish, simplifying the perceptual error to ( ( )) = ∏︁ * ( ( ))∏︁ 2 2 . This enables the optimization to consider swaps between any two pixels in the error image ( ). That image can be quickly rendered in advance, by invoking the render-function surrogate ′ with the input sample-set image.  take a similar approach, with swapping based on simulated annealing. Their empirically motivated optimization energy uses an explicit (Gaussian) kernel, but instead of computing an error image through a rendering surrogate, it postulates that the distance between two sample sets is representative of the difference between their corresponding pixel errors. Such a smoothness assumption holds for bi-Lipschitz-continuous functions. Their energy can thus be understood to compactly encode properties of a class of rendering functions.  adopt the approach of , but their energy function replaces the distance between sample sets by the difference between their corresponding pixel errors. The errors are computed using an explicit render-function surrogate. They optimize for a large number of simple surrogates simultaneously, and leverage a compact representation of Sobol sequences to also support progressive sampling. We relate these two prior works to ours more formally in supplemental Section 6, also showing how our perceptual error formulation can be incorporated into the method of .
The approach of Ahmed and Wonka [2020] performs on-thefly scrambling of a Sobol sequence applied to the entire image. Image pixels are visited in Morton Z-order modified to breaks its regularity. The resulting sampler diffuses Monte Carlo error over hierarchically nested blocks of pixels giving a perceptually pleasing error distribution. However, the algorithmic nature of this approach introduces more implicit assumptions than prior works, making it difficult to analyze.
Our theoretical formulation and optimization methods enable the construction of a-priori sampling masks in a principled way. For horizontal optimization, we recommend using our iterative algorithm (Alg. 2) which can bring significant performance improvement over simulated annealing; such speed-up was reported by  for dither-mask construction. Vertical optimization is an interesting alternative, where for each pixel one of several sample sets would be chosen; this would allow for varying the sample count per pixel. Note that the ranking-key optimization for progressive sampling of  is a form of vertical optimization.

Discussion
Our formulation expresses a-priori and a-posteriori optimization under a common framework that unifies existing methods. These two approaches come with different trade-offs. A-posteriori optimization utilizes sampled image information, and in a vertical setting requires no assumptions except for what is necessary for surrogate construction. It thus has potential to achieve high output fidelity, especially on scenes with complex lighting as it is oblivious to the shape and dimensionality of the rendering function, as first demonstrated by . However, it requires pre-sampling (also post-sampling in the horizontal setting), and the optimization is sensitive to the surrogate quality.
Making aggressive assumptions allows a-priori optimization to be performed offline once and the produced samples to be subsequently used to render any image. This approach resembles classical sample stratification where the goal is also to optimize sample distributions w.r.t. some smoothness expectations. In fact, our a-priori formulation subsumes the per-pixel stratification problem, since the perceptual error is minimized when the error image ( ) has both low magnitude and visually pleasing distribution. Two main factors limit the potential of a-priori optimization, especially on scenes with non-uniform multi-bounce lighting. One is the general difficulty of optimizing sample distributions in high-dimensional spaces. The other is that in such scenes the complex shape of the rendering function, both in screen and sample space, can easily break smoothness assumptions and cause high (perceptual) error.
To test the capabilities of our formulation, in the following we focus on the a-posteriori approach. In the supplemental document we explore a-priori optimization based on our framework. The two approaches can also be combined, e.g., by seeding a-posteriori optimization with a-priori-optimized samples whose good initial guess can improve the quality and convergence speed.

EXTENSIONS
Our perceptual error formulation (4) approximates the effect of the HVS PSF through kernel convolution. In this section we analyze the relationship between kernel and viewing distance, as well as the impact of the kernel shape on the error distribution. We also present extensions that account for the HVS non-linearities in handling tone and color.
Kernels and viewing distance. As discussed in Section 3.1, the PSF is usually modelled over a range of viewing distances. The effect of the PSF depends on the frequencies of the viewed signal and the distance from which it is viewed. Pappas and Neuhoff [1999] have found that the Gaussian is a good approximation to the PSF in the context of halftoning. They derived its standard deviation in terms of the minimum viewing distance for a given screen resolution: Here is the visual angle between the centers of two neighboring pixels (in degrees) for screen resolution (in 1⇑inches) and viewing distance (in inches). The minimum viewing distance for a given standard deviation and resolution can be contained via the inverse formula: = 2 tan 180 0.00954 2 −1 . Larger values correspond to larger observer distances; we demonstrate the effect of that in Fig. 3 where the images become increasingly blurrier. In Fig. 4a, we compare that Gaussian kernel to two well-established PSF models from the halftoning literature [Näsänen 1984;González et al. 2006]. We have found the differences between all three to be negligible; we use the cheaper to evaluate Gaussian in all our experiments.
Decoupling the viewing distances. Being based on the perceptual models of the HVS [Sullivan et al. 1991;, our formulation (4) assumes that the estimate and the reference are viewed from the same (range of) distance(s). The two distances can be decoupled by applying different kernels to and : = ∏︁ * − ℎ ℎ ℎ * ∏︁ 2 2 .
Minimizing this error makes appear from some distance similar to seen from a different distance ℎ ℎ ℎ . The special case of using a Kronecker delta kernel ℎ ℎ ℎ = , i.e., with the reference seen from up close, yields = ∏︁ * − ∏︁ 2 2 . This has been shown to have an edge enhancing effect [Anastassiou 1989;Pappas and Neuhoff 1999] which we show in Fig. 4b. We use ℎ ℎ ℎ = in all our experiments.
Tone mapping. Considering that the optimized image will be viewed on media with limited dynamic range (e.g., screen or paper), we can incorporate a tone-mapping operator into the perceptual error (4): Doing this also bounds the per-pixel error = ( ) − ( ), suppressing outliers and making the optimization more robust in scenes with high dynamic range. We illustrate this improvement in Fig. 4c, where an ACES [Arrighetti 2017] tone-mapping operator is applied to the optimized image. Optimizing w.r.t. the original perceptual error (4) yields a noisy and overly dark image compared to the tone-mapped ground truth. Accounting for tone mapping in the optimization through Eq. (13) yields a more faithful result.   Color handling. While the HVS reacts more strongly to luminance than color, ignoring chromaticity entirely (e.g., by computing the error image from per-pixel luminances) can have a negative effect on the distribution of color noise in the image. To that end, we can penalize the perceptual error of each color channel ∈ separately: where is a per-channel weight. In our experiments, we use an RGB space = {r, g, b}, set = 1, and use the same kernel = for every channel. Figure 4d shows the improvement in color noise over using greyscale perceptual error. Depending on the color space, the per-channel kernels may differ (e.g., YCbCr) [Sullivan et al. 1991].
As an alternative, one could decouple the channels from the input estimates and optimize each channel separately, assembling the results into a color image. In a vertical setting, this decoupling extends the optimization search space size from to ⋃︀ ⋃︀ .
Kernel shape impact. To test the robustness of our framework, we analyze kernels with spectral characteristics other than isotropic blue-noise in Fig. 5. We run our iterative pixel-swapping algorithm (Alg. 2) to optimize the shape of a white-noise input, which produces a spectral distribution inverse to that of the target kernel. The rightmost image in the figure shows the result of using a spatially varying kernel that is a convex combination between a low-pass Gaussian and a high-pass anisotropic kernel, with the interpolation parameter varying horizontally across the image. Our algorithm can adapt the noise shape well.

RESULTS
We now present empirical validation of our error optimization framework in the a-posteriori setting described in Section 5.2. We render static images and animations of several scenes, comparing our algorithms to those of .

Setup
Perceptual error model. We build a perceptual model by combining all extensions from Section 6. Our estimate-image kernel is a binomial approximation of a Gaussian [Lindeberg 1990]. For performance reasons and to allow smaller viewing distances we use a 3×3-pixel kernel with standard deviation = ⌈︂ 2⇑ (see Fig. 4a). Plugging this value into the inverse of Eq. (11), the corresponding minimum viewing distance is = 4792⇑ inches for a screen resolution of dpi (e.g., 16 inches at 300 dpi). We recommend viewing from a larger distance, to reduce the effect of our 3×3 kernel discretization. We use a Dirac reference-image kernel: ℎ ℎ ℎ = , and incorporate a simple tone-mapping operator that clamps pixel values to (︀0, 1⌋︀. The final error model reads: where ′ ′ ′ is the surrogate image whose construction we describe below. For dithering we convert RGB colors to luminance, which reduces the number of components in the error (15) to one.
Methods. We compare our four methods from Algs. 1 and 2 to the histogram and permutation of . For our vertical and horizontal iterative minimizations we set the maximum iteration count to 100 and 10 respectively. For error diffusion we use the kernel of Floyd and Steinberg [1976] and for dithering we use a void-and-cluster mask [Ulichney 1993]. For our horizontal iterative minimization we use a search radius = 1 and allow pixels to travel within a disk of radius = 1 from their original location in the dissimilarity metric. For the permutation method of  we obtained best results with tile size 8×8. (Our = 1 approximately corresponds to their tile size 3×3.) Rendering. All scenes were rendered with PBRT [Pharr et al. 2016] using unidirectional or bidirectional path tracing. None of the methods depend on the sampling dimensionality, though we set the maximum path depth to 5 for all scenes to maintain reasonable rendering times. The ground-truth images have been generated using a Sobol sampler with at least 1024 samples per pixel (spp); for all test renders we use a random sampler. To facilitate numerical-error comparisons between the different methods, we trace the primary rays through the pixel centers.
Surrogate construction. To build a surrogate image for our methods, we filter the per-pixel averaged input estimates using Intel Open Image Denoise [Intel 2018] which also leverages surface-normal and albedo buffers, taking about 0.5 sec for a 512 × 512 image. Recall that the methods of  utilize implicit surrogates.
Image-quality metrics. We evaluate the quality of some of our results using the HDR-VDP-2 perceptual metric [Mantiuk et al. 2011], with parameters matching our binomial kernel. We compute errordetection probability maps which indicate the likelihood for a human observer to notice a difference from the ground truth.
Additionally, we analyze the local blue-noise quality of the error image = ( ) − ( ). We split the image into tiles of 32 × 32 pixels and compute the Fourier power spectrum of each tile. For visualization purposes, we apply a standard logarithmic transform ln(1 + ⋃︀ˆ⋃︀) to every resulting pixel valueˆand compute the normalization factor per tile so that the maximum final RGB value within the tile is (1, 1, 1). Note that the error image is computed w.r.t. the ground truth and not the surrogate, which quantifies the blue-noise distribution objectively. The supplemental material contains images of the tiled power spectra for all experiments.
We compare images quantitatively via traditional MSE as well as a metric derived from our perceptual error formulation. Our perceptual MSE (pMSE) evaluates the error (15) of an estimate image w.r.t. the ground truth, normalized by the number of pixels and channels : pMSE = ⋅ . It generalizes the MSE with a perceptual, i.e., non-delta, kernel . Table 1 summarizes the results.

Rendering comparisons
All methods. Figure 6 shows an equal-sample comparison of all methods. Vertical methods select one of the 4 input samples per pixel; horizontal methods are fed a 2-spp average for every pixel, and another 2 spp are used to render the final image after optimization. Our methods consistently perform best visually, with the vertical iterative minimization achieving the lowest perceptual error, as corroborated by the HDR-VDP-2 detection maps. Error diffusion is not far behind in quality and is the fastest of all methods along with dithering. The latter is similar to Heitz and Belcour's histogram method but yields a notably better result thanks to using a superior surrogate and performing the thresholding as in the classical halftoning setting (see Section 5.2). Horizontal methods perform worse due to noisier input data (half spp) and worse surrogates derived from it, and also mispredictions (which necessitate re-rendering). Ours still uses a better surrogate than Heitz and Belcour's permutation and is also able to better fit to it. Notice the low fidelity of the 4-spp average image compared to our vertical methods', even though the latter retain only one of the four input samples for every pixel.
Vertical methods. In Fig. 7 we compare our vertical iterative minimization to the histogram sampling of . Both select one of several input samples (i.e., estimates) for each pixel. Our method produces a notably better result even when given 16× fewer samples to choose from. The perceptual error of histogram sampling does not vanish with increasing sample count. It dithers pixel intensity rather than pixel error, thus more samples help improve the intensity distribution but not the error magnitude. Figure 1 shows our most capable method: vertical iterative minimization with search space extended to the power set of the input samples (with size 2 4 − 1 = 15 for 4 input spp; see Section 4.3). We compare surrogate-driven optimization against the best-case resultoptimization w.r.t. the ground truth. Both achieve high fidelity, with little difference between them and with pronounced local blue-noise error distribution corroborated by the tiled power spectra.
Horizontal methods & animation. For rendering static images, horizontal methods are at a disadvantage compared to vertical ones due to the required post-optimization re-rendering. As  note, in an animation setting this sampling overhead can be mitigated by reusing the result of one frame as the initial estimate for the next. In Fig. 8 we compare their permutation method to our horizontal iterative minimization. For theirs we shift a voidand-cluster mask in screen space per frame and apply retargeting, and for ours we traverse the image pixels in different random order. We intentionally keep the scenes static to test the methods' best-case abilities to improve the error distribution over frames.
Starting from a random initial estimate, our method can benefit from a progressively improving surrogate that helps fine-tune the error distribution via localized pixel swaps. The permutation method operates in greyscale within static non-overlapping tiles. This prevents it from making significant progress after the first frame. While mispredictions cause local deviations from blue noise in both results, these are stronger in the permutation method's. This is evident when comparing the corresponding prediction images-the results of optimization right before re-rendering. The permutation's retargeting pass breaks the blocky image structure caused by tile-based optimization but increases the number of mispredictions.
The supplemental video shows animations with all methods, where vertical ones are fed a new random estimate per frame. Even without accumulating information over time, these consistently beat the two horizontal methods. The latter suffer from mispredictions under fast motion and perform similarly to one another, though ours remains superior in the presence of temporal smoothness. Mispredictions could be eliminated by optimizing frames independently and splitting the sampling budget into optimization and re-rendering halves (as in Fig. 6), though at the cost of reduced sampling quality.  Bottom row shows HDR-VPD-2 error-detection maps (blue is better, i.e., lower detection probability). The baseline 1-spp and 4-spp images exhibit large perceptual error, while our vertical iterative minimization achieves highest fidelity. Error diffusion produces similar quality. Dithering is as fast but shows smaller improvement over the baselines, yet significantly outperforms the similar histogram method. Our horizontal iterative optimization does better than the permutation method. Our methods also reduce MSE compared to the 4-spp baseline, even though they do not focus solely on per-pixel error (see Table 1). Additional comparisons. Figure 9 shows additional results from our horizontal and vertical minimization and error diffusion. We compare these to the permutation method of  which we found to perform better than their histogram approach on static scenes at equal sampling rates. For the horizontal methods we show the results after 16 iterations. Our methods again yield lower error, subjectively and numerically (see Tables 1 and 2).

DISCUSSION 8.1 Bias towards surrogate
While ultimately we want to optimize w.r.t. the ground-truth image, in practice we have to rely on a surrogate. In our experiments we use reasonably high-quality surrogates, shown in Fig. 12, to best demonstrate the capabilities of our framework. But when using a surrogate of low quality, fitting too closely to it can produce an estimate with artifacts. In such cases less aggressive fitting may yield lower perceptual error. To explore the trade-off, in Appendix B we augment the perceptual error with a term that penalizes deviations from the initial estimate init (which case of vertical optimization is obtained by averaging the input per-pixel estimates): The parameter ∈ (︀0, 1⌋︀ encodes our confidence in the surrogate quality. Setting = 1 reverts to the original formulation (15), while optimizing with = 0 yields the initial image estimate init . Optimizing w.r.t. this energy can also be interpreted as projecting the surrogate onto the space of Monte Carlo estimates in Ω Ω Ω, with control over the fitting power of the projection via . In Fig. 10, we plug the extended error formulation (16) into our vertical iterative minimization. The results indicate that the visually best result is achieved for different values of depending on the surrogate quality. Specifically, when optimizing w.r.t. the ground truth, the fitting should be most aggressive: = 1. Conversely, if the surrogate contains structural artifacts, the optimization should be made less biased to it, e.g., by setting = 0.5. Other ways to control this bias are using a more restricted search space (e.g., nonpower-set) and capping the number of minimization iterations of our methods. Note that the methods of  rely on implicit surrogates and energies and thus provide no control over this trade-off. We have found that their permutation method generally avoids tiling artifacts induced by their piecewise constant surrogate due to the retargeting step blurring the prediction image (shown in Fig. 8 zoom-ins); however, this blurring adds mispredictions which deteriorate the final image quality. Our methods provide better fits, target the error explicitly, and are much superior when the surrogate is good. With a bad surrogate, ours can be controlled to never do worse than theirs.

Denoising
Our images are optimized for eliminating error and preserving features when blurred with a given kernel. This blurring can be seen as a simple form of denoising, and it is reasonable to expect that the images are also better suited for general-purpose denoising  For high-quality surrogates (left and middle columns), the best result is achieved for values of close to 1. In contrast, strong structural artifacts (right column) call for lowering to avoid fitting too closely to the surrogate. The (subjectively) best image in each column is outlined in red.
than traditional white-noise renderings are Belcour and Heitz 2021]. However, we have found that obtaining such benefit is not straightforward.
In Fig. 11 we run Intel Open Image Denoise on the results from our vertical iterative minimization. On the left scene, the input samples ➀ have white-noise image distribution with large magnitude; feeding their per-pixel averages to the denoiser, it cannot reliably separate the signal from the noise and produces conspicuous artifacts. Using this denoised image ➁ as a surrogate for our optimization yields a "regularized" version ➂ of the input that is easier for the denoiser to subsequently filter ➃ . This process can be seen as projecting the initial denoised image back onto the space of exact per-pixel estimates (while minimizing the pMSE) whose subsequent denoising avoids artifacts. Note that obtaining this improved result requires no additional pixel sampling.
On the right scene in Fig. 11, the moderate input-noise level is easy for the denoiser to clean while preserving the faint shadow on the wall. Our optimization subsequently produces an excellent result which yields a high-fidelity image when convolved with the optimization kernel . Yet that same result is ruined by the denoiser which eradicates the shadow, even though subjectively its signalto-noise ratio is higher than that of the input image. Overall, the denoiser blurs our result ➂ aggressively on both scenes, eliminating not only the high-frequency noise but also lower-frequency signal not present in auxiliary input feature buffers (depth, normals, etc).
It should not be too surprising that an image optimized for one smoothing kernel does not always yield good results when filtered with other kernels. As an example, Fig. 5 shows clearly that the optimal noise distribution varies significantly across different kernels. While our kernel has narrow support and fixed shape, denoising kernels vary wildly over the image and are inferred from the input in order to preserve features. Importantly, modern kernel-inference models (like in the used denoiser) are designed (or trained) to expect mutually uncorrelated pixel estimates [Intel 2018]. This white-noiseerror assumption can also yield wide smoothing kernels that are unnecessarily aggressive for blue-noise distributions; we suspect  this is what hinders the denoiser from detecting features present in our optimized results whose pixels are highly correlated.
Our firm belief is that denoising could consistently benefit from error optimization, though that would require better coordination between the two. One avenue for future work would be to tailor the optimization to the kernels employed by a target denoiser. Conversely, denoising could be adapted to ingest correlated pixel estimates with high-frequency error distribution; this would enable the use of less aggressive smoothing kernels (see Fig. 3) and facilitate feature preservation. As a more immediate treatment, image features could be enhanced before or after our optimization to mitigate the risk of them being eliminated by denoising.

Performance and utility
Throughout our experiments, we have found that the tested algorithms rank in the following order in terms of increasing ability to minimize perceptual error on static scenes at equal sampling cost: histogram sampling, our dithering, permutation, our error diffusion, our horizontal iterative, our vertical iterative. The three lowestranked methods employ some form of dithering which by design assumes (a) constant image signal and (b) equi-spaced quantization levels shared by all pixels. The latter assumption is severely broken in the rendering setting where the "quantization levels" arise from (random) pixel estimation. Our vertical methods (dithering, error diffusion, iterative) are more practical than the histogram sampling of  as they can achieve high fidelity with a much lower input-sample count. Horizontal algorithms are harder to control due to their mispredictions which are further exacerbated when reusing estimates across frames in dynamic scenes.
Our iterative minimizations can best adapt to the input and also directly benefit from the extensions in Section 6 (unlike all others). However, they are also the slowest, as evident in Table 2 Practical utility. Our methods can enhance the perceptual fidelity of static and dynamic renderings as demonstrated by our experiments. For best results and maximum flexibility, we suggest using our vertical iterative optimization, optionally with the efficiency improvements mentioned above. Figure 10 illustrates that in practical scenarios (middle and right columns) this method can improve upon both the surrogate (top row) and the input-estimate average (bottom row) for a suitable value of the confidence parameter . For maximum efficiency we recommend using our vertical error diffusion. To obtain a surrogate, we recommend regularizing the input estimates via fast denoising or more basic bilateral or nonlocal-means filtering. Our optimization can then be interpreted as reducing bias or artifacts in such denoised images (see Fig. 10). Simple denoising of the result may yield better quality than traditional aggressive denoising of the input samples.
Progressive rendering. Our optimization methods produce biased pixel estimates through manipulating the input samples; this is true even for a-priori methods where the sampling is completely deterministic. Nevertheless, consistency can be achieved through a simple progressive-rendering scheme: For each pixel, newly generated samples are cumulatively averaged into a fixed set of per-pixel estimates that are periodically passed to the optimization to obtain an updated image. Each individual estimate will converge to the true pixel value, thus the optimized image will also approach the ground truth-with bounded memory footprint. Interestingly, convergence is guaranteed regardless of the optimization method and surrogate used, though better methods and surrogates will yield better starting points. Lastly, adaptive sampling is naturally supported by vertical methods as they are agnostic of differences in sample counts between pixels.

CONCLUSION
We devise a formal treatment of image-space error distribution in Monte Carlo rendering from both quantitative and perceptual aspects. Our formulation bridges the gap between halftoning and rendering by interpreting the error distribution problem as an extension of non-uniform multi-tone energy minimization halftoning.
To guide the distribution of rendering error, we employ a perceptual kernel-based model whose practical optimization can deliver improvements not achievable by prior methods given the same sampling data. Our model provides valuable insights as well as a framework to further study the problem and its solutions. A promising avenue for future research is to adapt even stronger perceptual error models. Prior work has already demonstrated a strong potential in reducing Monte Carlo noise visibility error using visual masking [Bolin and Meyer 1998;Ramasubramanian et al. 1999]. Robust metrics, other than squared ℒ 2 norm, can also be considered with possible nonlinear relationships.
Our framework could conceivably be extended beyond the human visual system, i.e., for optimizing the inputs to other types of image processing such as denoising. For such tasks, one could consider lifting the assumption of a fixed kernel to obtain an even more general problem where the kernel and sample distribution are optimized simultaneously (or alternatingly).
The first term in Eq. (17c) involves pixel permutations in the readily available estimated image ( ). In the second term we make an approximation that avoids rendering invocations: ∏︁Δ Δ Δ∏︁ 2 2 ≈ ∑ ( , ( )), where ( , ) measures the dissimilarity between the light-transport integrals of pixels and . Dropping the constant 2, we take the resulting bound as the optimization energy in Eq. (10b).

B SURROGATE CONFIDENCE CONTROL
Here we extend our perceptual error formulation to account for deviations of the surrogate image ′ ′ ′ from the ground truth . We introduce a parameter ∈ (︀0, 1⌋︀ that encodes our confidence in the quality of the surrogate and instructs the optimization how closely to fit to it. Given an initial image estimate init (the per-pixel estimate average for vertical optimization), we look to optimize for . We begin with an artificial expansion: ⌋︂ = ∏︁ * − ℎ ℎ ℎ * ∏︁ 2 (18a) = ∏︁(1 − )( * − * init ) + ( * − ℎ ℎ ℎ * ′ ′ ′ ) + (1 − )( * init − ℎ ℎ ℎ * ) + (ℎ ℎ ℎ * ′ ′ ′ − ℎ ℎ ℎ * )∏︁ 2 . (18b) Using the triangle inequality we then obtain the following bound: The second term on the right-hand side can be dropped as it is independent of the optimization variable . We bound the square of the first term using Jensen's and Young's convolution inequalities: We take this bound to be our optimization energy in Eq. (16), noting that the squared norm in the second term is the original energy with the surrogate ′ ′ ′ substituted for the ground truth . If a confidence map is available (e.g., as a byproduct of denoising), the minimization can be done with per-pixel control: In this supplemental document we discuss various details related to our general formulation from the main paper. We start with a description of the extension of our framework to the a-priori setting (Section 1). In Section 2 we describe a way in which textures can be accounted for in our horizontal approach, so that mispredictions due to multiplicative (and additive) factors are eliminated. In Section 3 we describe ways in which the runtime of iterative energy minimization methods can be improved considerably. Notably, an expression is derived allowing the energy difference due to trial swaps to be evaluated in constant time (no scaling with image size or kernel size).
In the remaining sections we analyze how current a-posteriori [Heitz and Belcour 2019] (Section 5) and a-priori ] (Section 6) state of the art approaches can be related to our framework. Interpretations are discussed, major sources of error are identified, and the assumptions of the algorithms are made explicit.

A-PRIORI OPTIMIZATION
We extend our theory to the a-priori setting and discuss the main factors affecting the quality. The quality of a-priori approaches is determined mainly by three factors: the energy, the search space, and the optimization strategy. We discuss each of those briefly in the following paragraphs.
Our energy. We extend the a-posteriori energy from the main paper in order to handle multiple estimators involving different integrands: 1 , . . . , , with associated weights 1 , . . . , : (1) In the above would typically be a low-pass kernel (e.g., Gaussian), and is the integral of the function used in the estimator . Through this energy a whole set of functions can be optimized for, in order for the sequence to be more robust to different scenes and estimators, that do not fit any of the considered integrands exactly. We note that the derived optimization in Section 3 below is also applicable to the minimization of the proposed energy.
Search space. The search space plays an important role for the qualities which the optimized sequences exhibit. A more restricted search space provides more robustness and may help avoid overfitting to the considered set of integrands.
For instance, sample sets may be generated randomly within each pixel. Then, their assignment to pixels may be optimized over the space of all possible permutations. This is the setting of horizontal methods. If additionally this assignment is done within each dimension separately it allows for an even better fit to the integrands in the energy (but may degrade the general integration properties of the sequence). The scrambling keys' search space in ] is a special case of the latter applied for the Sobol sequence.
Constraining the search space to points generated from lowdiscrepancy sequences provides further robustness and guarantees desirable integration properties of the considered sequences. Similarly to , we can consider a search space of Sobol scrambling keys in order for the optimized sequence to have a low discrepancy.
Ideally, such integration properties should arise directly from the energy. However, in practice the scene integrand cannot be expected to exactly match the set of considered integrands, thus extra robustness is gained through the restriction. Additionally, optimizing for many dimensions at the same time is costly as noted in ], thus imposing low-discrepancy properties also helps in that regard.
Finally, by imposing strict search space constraints a severe restriction on the distribution of the error is imposed. This can be alleviated by imposing the restrictions through soft penalty terms in the energy. This can allow for a trade-off between blue noise distribution and integration quality for example.
Progressive rendering. In order to make the sequence applicable to progressive rendering, subsets of samples should be considered in the optimization. Given a sample set for pixel we can decompose it in sample sets of 1, . . . , samples: ,1 ⊂ . . . ⊂ , ≡ . We denote the respective images of sample sets 1 , . . . , .
Then an energy that also optimizes for the distribution of the error at each sample count is: If , are set to zero for < then the original formulation is recovered. The more general formulation imposes additional constraints on the samples, thus the quality at the full sample count may be compromised if we also require a good quality at lower sample counts.
Choosing samples from for ,1 , . . . , , −1 (in each dimension) constitutes a vertical search space analogous to the one discussed in the main paper for a-posteriori methods. The ranking keys' optimization in ] is a special case of this search space for the Sobol sequence.
Adaptive sampling can be handled by allowing a varying number of samples per pixel, and a corresponding energy derived from the one above. Note that this poses further restrictions on the achievable distribution.
Optimization strategies. Typically the energies for a-priori methods have been optimized through simulated annealing . Metaheuristics can lead to very good minima especially if the runtime is not of great concern, which is the case since the sequences are precomputed. Nevertheless, the computation still needs to be tractable. The energies in previous works are generally not cheap to evaluate. On the other hand, our energies, especially if the optimizations in Section 3 are considered, can be evaluated very efficiently. This is beneficial for keeping the runtime of metaheuristics manageable, allowing for more complex search spaces to be considered. Implementation details. Implementation decisions for a renderer, such as how samples are consumed, or how those are mapped to the hemisphere and light sources, affect the estimator . This is important, especially when choosing for the described energies to optimize a sequence. It is possible that very small implementation changes make a previously ideal sequence useless for a specific renderer. It is important to keep this in mind when optimizing sequences by using the proposed energies and when those are used in a renderer.

TEXTURE DEMODULATION FOR HORIZONTAL OPTIMIZATION
Our iterative energy minimization algorithms (Alg. 1, Alg. 2, main paper) directly work with the original energy formulation, unlike error diffusion and dither matrix halftoning which only approximately minimize the energy. This allows textures to be handled more robustly compared to the permutation approach of Heitz and Belcour.
Comparison to demodulation. In the method of Heitz and Belcour the permutation is applied on the albedo demodulated image. This preserves the property that the global minimum of the implicit energy can be found through sorting. Translated to our framework this can be formulated as ( is a blue noise mask optimized for a kernel ): We have assumed that is zero, but we can also extend the method to handle an additive term as in our case. The more important distinction is that while the albedo demodulated image ′ is used in the permutation, it is never re-modulated ( ⊙ ⋅ is missing). Thus, this does not allow for proper handling of textures, even if it allows for modest improvements in practice. An example of a fail case consists of an image that is close to white noise. Then the error distribution will also be close to white noise due to the missing ⊙ ⋅ factor. More precisely, even if ( ′ ( )) − ′ ′ ′ is distributed as , this does not imply that ⊙ ( ′ ( )) − ′ ′ ′ will be distributed similarly. Dropping ⊙ ⋅ is, however, a reasonable option if one is restricted to sorting as an optimization strategy.
Sinusoidal textures. To demonstrate texture handling (multiplicative term ), in the top row of Fig. 1, a white-noise texture is multiplied with a sine-wave input signal: ( , ) = 0.5 * (1.0 + sin( + )) * ( , our horizontal iterative minimization algorithm can incorporate the albedo explicitly using the discussed energy. The bottom row demonstrates the effect of a non-flat signal on the error distribution (additive term ). Here is added to a sinewave input signal: ( , ) = 0.5 * (1.0 + sin( + )) + ( , ). The reference image is 0.5 * (1 + sin( + )). Our optimization is closer to the reference suggesting that our method can greatly outperform the current state of the art by properly accounting for auxiliary information, especially in regions with high-frequency textures.
Dimensional decomposition. The additive factor can be used to motivate splitting the optimization over several dimensions, since the Liouville-Neumann expansion of the rendering equation is additive ]. If some dimensions are smooth (e.g., lower dimensions), then a screen space local integrand similarity assumption can be encoded in (⋅, ⋅) and it will approximate Δ Δ Δ better for smoother dimensions. If the optimization is applied over all dimensions at the same time, this may result in many mispredictions due to the assumption being violated for dimensions in which the integrand is less smooth in screen space (e.g., higher dimensions). We propose splitting the optimization problem starting from lower dimensions and sequentially optimizing higher dimensions while encoding a local smoothness (in screen space) assumption on the integrand in (⋅, ⋅) (e.g., swaps limited to a small neighborhood around the pixel). This requires solving several optimization problems, but potentially reduces the amount of mispredictions. Note that it does not require more rendering operations than usual.

IMPROVING ITERATIVE-OPTIMIZATION PERFORMANCE
The main cost of iterative minimization methods is computing the energy for each trial swap, more specifically the required convolution and the subsequent norm computation. In the work of Analoui and Allebach an optimization has been proposed to efficiently evaluate such trial swaps, without recomputing a convolution or norm at each step, yielding a speed up of more than 10 times. The optimization relies on the assumption that the kernel is the same in screen space (the above optimization is not applicable for spatially varying kernels). We extend the described optimization to a more general case, also including spatially varying kernels. We also note some details not mentioned in the original paper.

Horizontal swaps
We will assume the most general case: instead of just swapping pixels, we consider swapping sample sets from which values are generated through . It subsumes both swapping pixel values and swapping pixel values in the presence of a multiplicative factor .
Single swap. The main goal is to evaluate the change of the energy due to a swap between the sample sets of some pixels , . More precisely, if the original sample set image is then the new sample set image is ′ such that ′ = , ′ = , and ′ = everywhere else. This corresponds to images: = ( ) and ′ = ( ′ ). The two images differ only in the pixels with indices and . Let: We will also denote the convolved images˜= * and˜′ = * ′ , and also =˜− . Specifically: We want to be able to efficiently evaluate = ∏︁˜′ − ∏︁ 2 − ∏︁˜− ∏︁ 2 , since in the iterative minimization algorithms the candidate with the minimum is kept. Using the above expressions for˜′ we rewrite as: where ,ℎ ( ) = ∑ ∈Z 2 ( − )ℎ( ) is the cross-correlation of and ℎ. We have reduced the computation of to the sum of only 4 terms. Assuming that , is known (it can be precomputed once for a known kernel) and that , is known (it can be recomputed after a sufficient amount of swaps have been accepted), then evaluating a trial swap takes constant time (it does not scale in the size of the image or the size of the kernel).
Multiple accepted swaps. It may be desirable to avoid recomputing , even upon accepting a trial swap. For that purpose we extend the strategy from  for computing , , where is the error image after swaps have been accepted: This implies:˜=˜+ ∑ =1 ( − + − ), and consequently: This allows avoiding the recomputation of , after every accepted swap, and instead, the delta on the +1-st swap with trial differences , is: where , is computed from , and , as derived in Eq. (17). This computation scales only in the number of accepted swaps since the last recomputation of , . We also note that , ( − ) evaluates to zero if − is outside of the support of , . Additional optimizations have been devised due to this fact ].
The corresponding expression for , is:

Multiple simultaneous updates
If the search space is ignored and the formulation is analyzed in an abstract setting it becomes obvious that the vertical approach corresponds to an update of a single pixel, while the horizontal approach corresponds to an update of two pixels at the same time. This can be generalized further. Let different pixels be updated per trial, and let there be trials that have been accepted since , has been updated. Let the pixels to be updated in the current trial be: +1 1 , . . . , +1 , and the accepted update at step be at pixels: 1 , . . . , . Let 0 = be the original image. We define the sequence of images: ∶ = −1 , ⇑ ∈ { 1 , . . . , } and otherwise let be given. Let = − −1 . Using the above notation we arrive at an expression for , : The change in the energy due to the + 1-st update is:

Implementation details
Leaky energy. Similar to the original paper , in our extension was computed for a "leaky energy" which extended the support of the image by convolution. That is reflected in the fact that the sums are over Z 2 . To rectify this, the sum needs to be limited to the support of . This would require clamped sums of the cross-correlation to be evaluated, which can also be precomputed but requires extra memory. The same holds for the cross-correlation with , where clamped terms are required near the image boundary.
Reflecting boundary conditions. Another desirable property may be a convolution such that it acts on the image extended to be reflected at the boundaries -this avoids artifacts near the borders. This can be achieved by including the relevant terms including pixels for which the kernel is partially outside of the support of . Care must be taken when expressing˜, however, since it may include the same updated pixel numerous times (especially if it is near the border). The same ideas apply for a toroidally extended convolution.
Further optimizations. Various other strategies have been proposed in the literature for improving the runtime of iterative error minimization approaches for halftoning.
In our algorithms we usually use a randomized initial state, however, it is possibly to initialize the algorithms with the result of a dither matrix halftoning algorithm or error diffusion algorithm which would result in faster convergence .
Another strategy involves partitioning the image in blocks. Instead of updating the pixels in raster or serpentine order, the blocks are updated simultaneously by keeping only the best update per block in each iteration. This has been reported to run 10+ times faster [Lieberman and Allebach 1997]. In the same paper [Lieberman and Allebach 1997], approximating the kernel with box functions has been proposed, yielding a speed up of 6 times. Similarly, if the kernel is separable or can be approximated by a separable kernel, the convolution can also be made considerably faster. A speed-up of an additional 30 times has been reported in  through the usage of a GPU.
Finally, several heuristics related to the order in which pixels are iterated over have been proposed in [Bhatt et al. 2006].
Our goal is to evaluate the change in the energy due to the update: In the above , ( ) = ∑ ∈Z 2 , , may be precomputed for every , which yields a function with support supp( , ) = ⋃ supp( ), and , ( ) = ∑ ∈Z 2 , can also be recomputed after enough updates have been accepted.
Multiple accepted updates. Let a set of accepted updates results in the differences: { 1 , . . . , }. And let be the error image after the updates. We derive an expression for the efficient evaluation of , : An efficient computation of , can then be achieved if the function , ( , ) = ∑ ∈Z 2 , , is precomputed. Then, at step + 1 the change in energy is: Multiple simultaneous updates. We derive an expression where an update consists of changing pixels simultaneously, and we assume that such updates have been accepted previously. We denote the differences of the pixels in update : { 1 , . . . , }. The expression for the change in the energy is given as: Where , ( , ) = ∑ ∈Z 2 , , is assumed to be precomputed, and , can be computed as:

RELATIONSHIP TO PREVIOUS WORK
We show that the recent publications  on blue noise error distribution for path tracing, can be seen as special cases in our framework. This allows for a novel analysis and interpretation of the results in the aforementioned works. We also state the necessary assumptions and approximations necessary to get from our general formulation to the algorithms presented in the papers.
Classification. The proposed techniques can be divided into apriori ] and a-posteriori . The main difference is that for a-priori techniques broad assumptions are made on the integrand without relying on information from renderings of the current scene. The cited a-priori approaches describe ways for constructing offline optimized point sets/sequences. We denote the method in  as BNDS (blue-noise dithered sampling), the method in ] as HBS (Heitz-Belcour Sobol), and the histogram and permutation method in  as HBH and HBP respectively (Heitz-Belcour histogram/permutation).
Energy. HBH/HBP both rely on a blue noise dither matrix optimized while using a Gaussian kernel (through void-and-cluster [Ulichney 1993]). This kernel corresponds to the kernel in our framework . The optimization of this dither matrix happens offline unlike in our iterative energy minimization algorithms. This imposes multiple restrictions while allowing for a lower runtime. On the other hand, the dither matrices in HBS and BNDS are optimized with respect to empirically motivated energies that cannot be related directly to what is used as energy in HBH and HBP. In the case of BNDS the energy does not even introduce an implicit integrand, and instead it is devised to represent a whole class of integrands.
We propose to substitute those empirically motivated energies with a modified version of our energy. This allows an intuitive interpretation and relating a-posteriori approaches to a-priori approaches.
Search space. Another notable difference constitute the search spaces on which the different approaches operate. HBH selects a subset from a set of precomputed samples in each pixel, HBP permutes the assignment of sample sets to pixels, BNDS directly modifies the set of samples in each pixel, and HBS considers a search space made up of scrambling and ranking keys for a Sobol sequence. Working on the space of scrambling and ranking keys guarantees the preservation of the desirable integration qualities of the Sobol sequence used, and it should be clear that other methods can also be restricted to such a space. Clearly, a search space restriction diminishes the achievable blue noise quality. On the other hand, it makes sequences more robust to integrands for which those were not optimized.

A-POSTERIORI APPROACHES
In this section we analyze the permutation based approach (HBP) and the histogram sampling approach (HBH) proposed in . The two methods can be classified as dither matrix halftoning methods in our framework, that operate on a horizontal and vertical search space respectively. We make the approximations and assumptions necessary to get from our general formulation to HBP/HBH explicit.
We also note that a-posteriori methods lead to solutions that adapt to the current render by exploiting known information (e.g. previously rendered data, auxiliary buffers). They can generally produce better results than a-priori methods.
Both HBP and HBH rely on a blue noise dither matrix . Let be the optimized blue noise dither matrix resulting from the minimization of ( ) = ∏︁ * ∏︁ 2 2 over a suitable search space. The kernel is the one used to generate the blue noise images for HBP/HBH. That is, the Gaussian kernel in the void-and-cluster method [Ulichney 1993]. Our analysis does not rely on the kernel being a Gaussian, or on the void-and-cluster optimization, this is simply the setting of the HBP/HBH method. In the more general setting any kernel is admissible.

Sorting step for the permutation approach
The permutation approach  consists of two main parts: sorting (optimization), and retargeting (correcting for mispredictions). The sorting step in HBP can be interpreted as minimizing the energy: A global minimum of the above energy is achieved for a permutation that matches the order statistics of and . Thus our goal would be to get from the minimization of: to the minimization of Eq. (51) over a suitable search space (in practice it is limited to permutations within tiles). We successively bound the error, while introducing the assumptions implicit to the HBP method. The bounds are not tight, however, the different error terms that we consider illustrate the major sources of error due to the approximation of the more general energy (Eq. (52)) with a simpler one (Eq. (51)). The substitution of the kernel convolution * ⋅ by a difference with a blue noise mask restricts the many possible blue noise error distributions towards which ( ( )) can go with a single one: . A global minimizer of the new simplified energy can thus be found by just sorting.
The closer the distributions of ( ( )) and , > 0 are locally, the lower this restriction error can be made. Notably, for a close to linear relationship between the samples and the integrand, and sufficiently many pixels, ( ( )) and can be matched closely in practice. A different way to reduce the approximation error is to introduce a sufficient amount of different blue noise images and pick the one that minimizes the error. We start with the original energy (Eq. (52)) and bound it through terms that capture the main assumptions on which the model relies: In the above, 2 is taken over the space of all strictly monotonically increasing functions, and > 0 is a real value used to provide an amplitude matching between ( ( )) and (this allows for the second term to go to zero as the pointwise error goes to zero).

Third error term.
We note that is precomputed offline in order to approximately minimize ( ) = ∏︁ * ∏︁ 2 . Thus, the third term reflects the quality of the blue noise achieved with respect to in the offline minimization. This error can be made small without a performance penalty since the optimization is performed offline. We factor out a multiplicative scaling factor > 0 in the blue noise quality term, to allow for the second term to go to zero. With this change, we can consider to be normalized in the range (︀−1, 1⌋︀ and we can encode the scaling in . 5.1.2 Second error term. The second term reflects the error introduced by substituting a large search space (many local minima) with a small search space. It introduces the first implicit assumption of HBP by relating the first and third error terms (by using 2 and respectively) through the second error term. The assumption is that there exists a permutation for which ( ( )) can be made close to , which would make the second term small. This holds in practice if the pixel-wise error is zero on average (unbiased estimator within each pixel), and we have a sufficiently large resolution/tiles: which results in a higher probability that pixels from ( ( )) can match well. Then the term ∏︁ ∏︁ 1 ∏︁ 2 ( ) − ∏︁ 2 can be made small. We note that this is a generalization of the third optimality condition in  (correlation-preserving integrand) since an integrand linear in the samples can also better match provided enough pixels. For a linear integrand the optimal 2 is also a linear function (ideal correlation between samples and integrand). The main difference between a linear integrand and a nonlinear/discontinuous one, is the amount of sample sets/pixels necessary to match 2 ( ) well, given an initial white noise samples' distribution. So in practice there are 4 factors directly affecting the magnitude of the second term: the number of considered blue noise images, the size of the tiles, the correlation between samples and integrand (accounted for by 2 ), the bias/consistency of the estimators.
We note that the number of considered pixels depends on the tile size in HBP, and the practical significance of this has been demonstrated through a canonical experiment in the main paper.

First error term.
Before we proceed we need to further bound the first error term by substituting ( ( )) by ( ( )). As discussed in the main paper, this is achieved by introducing a difference term Δ Δ Δ( ) = ( ( )) − ( ( )), and then ⌋︂ is recovered. The error there can be made arbitrarily small through 2 (it is accounted for in the second term). Thus we only need to study the remaining error due to Δ Δ Δ. In the case of HBP, Δ Δ Δ is approximated by non-overlapping characteristic functions in each tile ( ( , ) = ∞, for , in different tiles). This means that the approximation error is zero within each tile if the integrands are the same within the tile and permutations act only within the tile, since Δ Δ Δ( ) = 0 0 0. On the other hand, if this assumption is violated, mispredictions occur, usually resulting in white noise.
5.1.4 Δ Δ Δ term. HBP partitions screen space into a several tiles ℛ 1 , . . ., ℛ , and permutations are only over the pixel values in a tile. Having the partition induced by the tiling we can bound the first term: Since additionally the permutations are optimized for the pixel values instead of the sample sets (which saves re-rendering operations), then there is an assumption that within each tile ℛ the following holds (we denote = ⋃︀ ℛ ) ∶ ( ( )) = ( ( )).
Consequently it follows that = , ∀ , ∈ ℛ . This assumption can be identified with the 4-th optimality condition proposed in : screen-space coherence. As discussed, the search space restriction to the tiles corresponds to an approximation of the Δ Δ Δ term in our framework by characteristic functions: ( , ) = ∞, ∈ ℛ , ⇑ ∈ ℛ and ( , ) = 0, , ∈ ℛ . To account for the actual error when the assumption is violated we introduce an additional error term per tile: Δ Δ Δ = ( ( )) − ( ( )), then we have the bound: This means that even if all of the previous error terms are made small, including ∏︁ ( ( ))− 2 ( )∏︁ 2 , the error may still be large due to ∏︁Δ Δ Δ∏︁ 2 . We refer to a large error due to the delta term as misprediction -that is, a mismatch between the predicted error distribution from the minimization of ∏︁ ( ( )) − 2 ( )∏︁ 2 and the actual error distribution resulting from the above permutation applied to ( ( )). The best way to identify mispredictions is to compare the predicted image ( ( )) and the image rendered with the same permutation for the sample sets ( ( )). A misprediction occurring means that the assumption made to approximate Δ Δ Δ was incorrect ( Δ Δ Δ ≠ 0 0 0 for some tile ℛ ), equivalently the optimality condition of screen-space coherence is not satisfied.  Fig. 2. Here we showcase the effect of tile size on the quality of blue noise. We also demonstrate the effect of a larger search neighborhood in our optimization Alg. 2 (main paper). For our case, we consider disk neighborhoods so that they are contained within Heitz and Belcour's tiles in terms of size, but they can also overlap due to our formulation. From left-to-right, the input white noise texture is optimized using our relocation algorithm. The last two columns are from  method. The corresponding power spectra of these optimized images (128 × 128) are also shown.
Avoiding mispredictions. In practice mispredictions often occur for larger tile sizes, since it is hard to guarantee that the integrand remains similar over each tile. On the other hand, larger tiles allow for a better blue noise as long as Δ Δ Δ = 0 in each tile, thus larger tiles are desirable. The method fails even more often near edges, since even for small tile sizes it allows swapping pixels over an edge. A straightforward improvement involves partitioning the domain by respecting edges. More involved methods may take into account normals, depth, textures, etc.
We have dropped the term since it does not affect the set of minimizers ( is assumed constant in each tile). As discussed in Eq. (51), a global minimum is given by matching the order statistics of to the order statistics of 2 ( )) (we note that the order statistics of do not change from the application of 2 since it is a strictly increasing function). This is equivalent to performing the sorting pass described in . A minor optimization would be to pre-sort and instead store the sorted indices.
Tiling effect. In Fig. 2 we compare the effect of the tile size. In our approach, the "tiles" can be defined per pixel, can have arbitrary shapes, and are overlapping, the last being crucial for achieving a good blue noise distribution. We consider white-noise with mean 0.5 (which is an ideal scenario for Heitz and Belcour's method) and compare various tile sizes. For a fair comparison, our tile radius corresponds similar tile-size in the permutation [2019] approach.
The power-spectrum profiles confirm the better performance of our method. Retargeting [2019] cannot improve the quality of the permutation approach either, since no misprediction can occur (Δ Δ Δ = 0). The adverse effect of tiling is exacerbated in practice since, for images which are not smooth enough in screen space, tiles of smaller sizes need to be considered.
Custom surrogate. The term does not need to be assumed constant in fact. If it is assumed constant, that is equivalent to picking a tile-constant surrogate, however, a custom surrogate may be provided instead. Then one would simply minimize the energy: The energy has a different minimizer than the original HBP energy, but the global minimum can be found efficiently through sorting once again.

HBP retargeting
The retargeting pass in HBP achieves two things. It introduces new possible target solutions through new blue noise images, and it corrects for mispredictions. The first is not so much a result of the retargeting, as it is of varying the blue noise image every frame. Ideally several blue noise images would be considered in a single frame, and the best image would be chosen per tile (in that case one must make sure that there are no discontinuities between the blue noise images' tiles) in order to minimize the second term in Eq. (53). Instead, in HBP this is amortized over several frames.
The more important effect of retargeting is correcting for mispredictions, by transferring the recomputed correspondence between sample set and pixel value (achieved through rerendering) to the next frame. This allows reducing the error due to the approximation of Δ Δ Δ (when the piecewise-tile constancy assumption on the integrand is violated). Note however, that this is inappropriate if there is a large temporal discontinuity between the two frames. Implementation details. Retargeting requires a permutation that transforms the blue noise image in the current frame into the blue noise image of the next frame . This permutation is applied on the optimized seeds to transfer the learned correspondence between sample sets and pixel values to the next frame. Implicitly, this transforming permutation also relies on a screen space integrand similarity assumption, since there is no guarantee that the corresponding values from the swap will match, possibly incurring a misprediction once again (it can be modeled by an additional Δ Δ Δ term). In HBP ] the maximum radius of travel of each pixel in the permutation is set to 6 pixels. This has a direct effect on the approximation of Δ Δ Δ, as the travel distance of a pixel is allowed to extend beyond the original tile bounds. In the worst case scenario a pixel may allowed to travel a distance of ⌉︂ 2 + 2 + 6 pixels, where , are the dimensions of the tiles. An additional error is introduced since the retargeting pass does not produce the exact blue noise image used in the next frame, but some image that is close to it . This seems to be done purely from memory considerations since it allows one blue noise image to be reused by translating it toroidally each frame to produce the blue noise image for the next frame.
Relationship to our horizontal approach. Our horizontal approach does not require a retargeting pass. It can directly continue with the optimized sample sets and pixel values from last frame. There is also no additional travel distance for a matching permutation as in retargeting, which further minimizes the probability of misprediction. Thus, it inherently and automatically produces all of the advantages of retargeting while retaining none of its disadvantages.

Histogram sampling approach
The histogram sampling approach from Heitz and Belcour can be interpreted as both a dithering and a sampling method. We study the dithering aspect to better understand the quality of blue noise achievable by the method.
Algorithm analysis. The sampling of an estimate in each pixel by using the corresponding mask value to the pixel can be interpreted as performing a mapping of the mask's range and then quantizing to the closest estimate. In HBH each estimate is equally likely to be sampled (if a random mask is used), which implies a transformation that maps equal parts of the range to each estimate. Let ,1 , . . ., , be the greyscale estimates in pixel sorted in ascending order. Let the range of the blue noise mask be in [0,1]. Then the range is split into equal subintervals: (︀0, 1 ), . . . , (︀ −1 , 1⌋︀ which respectively map to (︀ 1 , 1 + 2 2 ), . . . , ⌋︀. If the quantization rounds to the closest estimate, then the above mapping guarantees the desired behavior. We note that since the estimates in each pixel can have different values, the mapping for each pixel may be different. We will denote the above mapping through . Then the mapping plus quantization problem in a pixel may be formulated as: Note that the minimization in each pixel is independent, and it aims to minimize the distance between the estimates and the remapped value from the blue noise mask. If the set of estimates are assumed to be the same across pixels, and are also assumed to be spaced regularly, then is only a linear remapping, which effectively transfers the spectral properties of onto the optimized image. Notably, the former is the screen-space coherence assumption from HBP, while the latter is the correlation-preserving integrand assumption. Thus we have seen that for optimal results the HBH method relies on exactly the same assumption as the HBP method (while our vertical iterative minimization approach lifts both assumptions).
Disadvantages. One of the key points is that the error distribution and not the signal itself ought to ideally be shaped as . This is actually the case even in the above energy. From the way was chosen it follows that the surrogate is equivalent to (0.5) which can be identified as the image made of the median of the sorted estimates within each pixel. This is the case since if the target surrogate of (during the offline optimization) was assumed to be 0.5, then after the mapping it is (0.5). Generally, this is a very bad surrogate in the context of rendering, and it generally increases the error compared to the averaged image, making the method impractical.
Another notable disadvantage is that all estimates are considered with an equal weight. This means that outliers are as likely to be picked as estimates closer to the surrogate. This results in fireflies appearing even when those were not present in the averaged imaged. Compared to classical halftoning, where only the closest lower and upper quantization levels are considered, HBH does not minimize the magnitude of the error to the surrogate.
Finally, the two assumptions of: screen-space coherence and correlation-preserving integrand, generally do not hold in practice. Estimates cannot be assumed to match between pixels (especially if samples are taken at random), and they cannot be assumed to be uniformly distributed, which implies that is not linear. This greatly impacts the quality of the result, especially if it is compared to adaptive approaches such as our vertical error diffusion approach and our iterative minimization techniques (see the experiments in the main paper).
Generalization. The method can be generalized to take a custom surrogate instead of the one constructed by the median of the estimates within each pixel. This is achieved by splitting the per pixel set of estimates into two parts: (greyscale) estimates greater than the value of the (greyscale) surrogate in the current pixel, and estimates lower than it. Then the mapping for the current pixel maps values in (︀0, 0.5) to the lower set, and values in (︀0.5, 1⌋︀ to the higher set, such that (0.5) = . The original method is recovered if the surrogate is chosen to be the implicit one for the original histogram sampling method and if the appropriate corresponding mapping is kept.
The approach can be extended further by setting different probabilities for the different estimates. The original histogram sampling method correspond to setting the same probability for sampling every estimate, equivalently: equal sized sub-intervals from (︀0, 1⌋︀ map to each estimate. Classical dither matrix halftoning can be interpreted as setting an equal probability for the closest to the surrogate upper and lower estimates, while every other estimate gets a zero probability. Equivalently: equal sub-intervals from (︀0, 1⌋︀ map to the two aforementioned estimates while no part of the interval maps to the remaining estimates. Generally a custom probability can be assigned to each estimate: 1 , ..., , by having the intervals (︀0, 1 ), ..., (︀∑ −1 =1 , 1⌋︀ map to 1 , ..., (after quantization). We note that an unbiased image can be recovered only if there is a map to every estimate.

A-PRIORI APPROACHES
We discuss current state of the art a-priori approaches  and their relation to our framework, as well as insights regarding those.

HBS
In Heitz et al.'s work, a scrambling energy and a ranking energy have been proposed (note that those energies are maximized and not minimized): The upper indices in 1 , 2 indicate that the two energies are evaluated with different subsets of the sample set in the pixel . The are taken from an arbitrary set of functions (in the original paper those are random Heaviside functions). The described form of the energies has been partially motivated by the energy in . This does not allow for a straightforward interpretation or a direct relation to the (implicit) energy used for a-posteriori approaches in ].
Scrambling energy. We modify in order to relate it to the energy in our framework and to provide a meaningful interpretation: We have relaxed the Gaussian kernel to an arbitrary kernel and absorbed it into the norm. More importantly we have removed the heuristic dependence of error terms on their neighbors, and instead the coupling happens through the kernel itself. Finally, we have introduced weights 1 , . . . , that allow assigning different importance to different integrands. Thus, this is a weighted average of our original energy applied to several different integrands, matching our a-priori approach (Eq. (1)). Through this formulation a direct relationship to the a-posteriori methods can be established, and it can be motivated in the context of both the human visual system and denoising. Particularly, the scrambling energy ′ is over the space of scrambling keys, which allow permuting the assignment of sample sets. This is in fact the horizontal setting from our formulation in the main paper. The space can be extended further if the scrambling keys in each dimension are different (as in HBS). The same can be done in a-posteriori methods, if the optimization is performed in each dimension as discussed in Section 2.
Ranking energy. The ranking keys in HBS describe the order in which samples are consumed. This is useful for constructing progressive a-priori methods. Notably, the order in which samples will be introduced can be optimized. Having a sequence of sample sets in each pixel: ,1 ⊂ . . . ⊂ , ≡ and respectively the images formed by those: 1 , . . . , , the progressive energy may be constructed as: The quality at a specific sample count corresponding to is controlled through the weight . The original energy maximizing the quality of the full set is retrieved for ( 1 , . . . , −1 , ) = (0, . . . , 0, 1). Since the sample sets , . . . , , are optimized by choosing samples from this can be seen as a vertical method. Finally, the ranking keys can also be defined per dimension, which can be related to a-posteriori methods through the suggested dimensional decomposition in Section 2.

Blue-noise dithered sampling energy
In Georgiev and Fajardo's work, in order to get an optimized (multichannel) blue noise mask, the following energy has been proposed: which bears some similarity to the weights of a bilateral filter. In the above , are pixel coordinates, and , are -dimensional vectors associated with , . Let the image formed by those vectors be . The energy aims to make samples , distant (∏︁ − ∏︁ must be large) if they are associated with pixels which are close (∏︁ − ∏︁ is small).
Relation to our framework. Even though the energy is heuristically motivated, we can very roughly relate it to our framework. The above energy implicitly assumes classes of integrands 1 , ..., , such that close samples , are mapped to close values , ( ), , ( ), and distant samples are mapped to distant values. Notably, the form of the energy does not change over screen-space, so the same can be implied about the integrands. One such class is the class of bi-Lipschitz functions. The bound can be used to relate a modified version of the original energy, to an energy of the form: (69) Thus, the original energy can indeed be interpreted as reasonable for a whole class of sufficiently smooth integrands, instead of an energy that works very well with one specific integrand.
A similar thing can be achieved in our framework, if the weighted energy is considered: The kernel can be a Gaussian with standard deviation , as in the original energy, or it can be relaxed to an arbitrary desired kernel. 1 , . . . , are representative integrands that satisfy the discussed smoothness requirements, and are associated weights assigning different importance to the integrands. Finally, the reference images are given by the integrals = ∫ (︀0,1⌋︀ ( ) .
It should be clear that this is a weighted average constructed from the standard energy in our framework applied to a set of integrands.
There are a number of benefits of such an explicit formulation. Most importantly, it allows for a-priori methods to be studied in the same framework as a-posteriori approaches. Additionally, explicit control is provided over the set of integrands and the kernel in a manner that allows for a straightforward interpretation.
Perceptual quality trade-off. While the energy of Georgiev and Fajardo is able to account for many different integrands, this is achieved at the cost of the perceptual quality of the produced patterns. We illustrate this in Fig. 4 by considering a constraint where 25% of all pixels have an error of +1 and 25% of all pixels have an error of -1.
For the experiment an initial white noise image is permuted using a brute force optimization with our energy from the main paper and the energy of Georgiev and Fajardo. One can see that the pattern resulting from our energy always decays faster under convolution. This can be explained by the fact that the bilateral filter-like energy forces nearby pixels to be as different as possible. This doesn't necessarily lead to the best results under convolution illustrated by Figure 5, but it is necessary in the setting of a-priori methods since not much information is assumed regarding the integrand.
We consider a more realistic example in Fig. 6 where the underlying signal is a sine function with vertically increasing frequency. We first degrade the signal with uniform white noise. To optimize the error distribution, we use our Kronecker kernel energy extension (eq. 11 from the main paper where ℎ = ) that is given by: where simply clamps values to (︀0, 1⌋︀. The result with our energy function matches better the original signal. This is perfectly in line with all of our results on realistic scenes presented in the main paper and the supplemental HTML.

Blue-noise dithered sampling algorithm
The second contribution of Georgiev and Fajardo's work is a sampler which relies on an image optimized with Eq. (68) and uses it to achieve a blue noise distribution of the rendering error. We summarize the algorithm and discuss some details related to it.
The sequence can be constructed by using various samplers (e.g., random, low-discrepancy, blue-noise, etc.). The construction of the new points for pixel can be interpreted either as toroidally shifting the sequence by or equivalently as toroidally shifting the sequence { , . . . , } by .
The sequences constructed within each pixel are used to estimate the integral in the usual manner. Since a finite number of dimensions are optimized the suggestion is to distribute the constructed sequences over smoother dimensions, while other dimensions may use a standard sampler.  . 3. A comparison illustrating that even a sampling sequence formed by a stack of blue noise images (Ours) yields a good distribution (note the tiled error spectra). The integration error is higher however, degrading the quality. This is the case because the assumed integrand is far from linear in each dimension (see Extension in Section 6.3). The images use 4 samples per pixel, and the degradation of the spectral properties with the number of samples is clear for  and even ], while it is not so much the case for Ours. This demonstrates that different methods offer a different trade-off between integration error and distribution for arbitrary integrands. Constraining the search space to using toroidal shifts or scrambling and ranking keys restricts the achievable blue noise distribution.

Error
Convolved error Absolute diff. Power spectrum  White noise Fig. 4. We show an example demonstrating how our energy (top row) forms clusters where required so that the convolved error (second column) produces the best cancellation effect. The first column shows error images.

Ours
Ours would converge to a grey (reference) image faster compared to the one using the energy in Eq. (68). The convolved images in the second column show the same behavior. The third column shows the absolute difference between the convolved error and the reference grey image (darker is better). The fourth column shows the error power spectra, with ours showing much better blue-noise characteristics than others.
Effect of the toroidal shift. Let us consider a linear one-dimensional integrand ( ) = + that does not vary in screen space, and a sequence with a single point . Furthermore, if we assume = 0, then the error is given by: Since does not vary in screen space, then also does not. Then the power spectrum of the error (excluding the DC) matches the Error Convolved error Absolute difference Ours  Fig. 5. We consider an example with 2 error pixels (+1 and -1). The first column shows the error images, the second column shows this error convolved with a gaussian kernel, and the third column shows the difference between the convolved error and the reference (constant) greyscale image. In the top row, our energy clusters these pixels such that they can cancel out each other's contribution under convolution. Georgiev and Fajardo's energy in the bottom row pushes these pixels farther away. The corresponding absolute difference (convolved error − constant grey image) images in the third column demonstrate that our energy makes the error converge faster to the constant greyscale image (darker is better).
power spectrum of up to the multiplicative factor 2 . Then, under the assumption that the integrand is linear, does not vary in screen space, and there is no toroidal shift, the power spectral properties of are transferred ideally to the error.
On the other hand, if is chosen to be non-zero, then the spectral characteristics of the image (( + ) mod 1) will be transferred instead. We have empirically verified that even with a very good quality blue noise image the toroidal shift degrades its quality due to the introduced discontinuities. Thus, even in the ideal case Initial Georgiev [2016] Ours Reference Fig. 6. A more realistic test with kernel using = 1⇑ ⌋︂ 2. The signal is a sine function that increases in frequency along the vertical axis. Our method handles tone mapping and preserves well both the lower and higher frequencies present in the signal.
of a constant in screen space linear 1-D integrand, toroidal shifts degrade the quality.
Effect of using multiple samples. Let us consider the same integrand ( ) = + , which we have identified as being ideal for transferring the spectral characteristics of to the error. And let us assume that we are given several samples: = { 1 , ..., }, and we have constructed the sample set image through toroidal shifts with . Then the error is: The power spectrum of the error thus matches the power spectrum of the image = ∑ =1 , (excluding the DC) up to a multiplicative factor. For a random point sequence the more points are considered, the closer to white noise becomes. This is further exacerbated by the discussed discontinuities introduced by the toroidal shifts.
Extension. We have argued that both toroidal shifts and increasing the number of samples has a negative effect on transferring the spectral properties of even in an ideal scenario. Naturally the question arises whether this can be improved. Our proposal is the direct optimization of point sets without the application of a toroidal shift.
For the discussed example this entails constructing a sequence of images 1 , . . . , such that = ∑ =1 is a blue noise image.