I describe a simple method for constructing a sequence of points, that is based on a low discrepancy quasirandom sequence but exhibits enhanced isotropic blue noise properties. This results in fast convergence rates with minimal aliasing artifacts.     Introduction Quasirandom low discrepancy sequences are useful to create distributions that appear less regular then lattices, but more regular than simple random distributions (see figure 2). They play a large role in many numerical computing applications, including physics, finance and in more recent decades computer graphics. The  classic paper on how using quasirandom sampling over uniform sampling can improve 2D images is “Stochastic Sampling in Computer Graphics” [Cook, 1986]. In more recent decades, quasirandom sampling has become a major tool in 3D rendering. Rendering (and samplers) can utilize either the CPU or GPU for calculations. For an excellent primer on their use in 3D rendering, see “Anti-aliased Low Discrepancy Samplers for Monte Carlo Estimators in Physically Based Rendering” [Perrier, 2108]. In this paper, she describes the task of 3D rendering (see figure 3) as follows: “When displaying a 3-D object on a computer screen, this 3-D scene is turned into a 2-D image, which is a set of organized colored pixels. We call Rendering the process that aims at finding the correct color to give those pixels. This is done by simulating and integrating all the light,….we numerically approximate the equation by taking random samples in the integration domain and approximating the integrand value using Monte Carlo methods….A stochastic sampler should minimize the variance in integration to converge to a correct approximation using as few samples as possible. Many different samplers exist that we roughly regroup into two families: –  Blue Noise samplers, that have a low integration variance while generating unstructured point sets. The issue with those samplers is that they are often slow to generate a point set – Low Discrepancy samplers, that minimize the variance in integration and are able to generate and enrich a point set very quickly. However, they present a lot of structural artifacts when used in Rendering.”  In recent years, efforts have been made to combine to the anti-aliasing advantages of jittered sampling, with the fast convergence rates (low variance) associated with low discrepancy quasirandom sequences [Christensen, 2018]. With this same goal, I present a new and different method to construct an infinite point sequence, that is extremely  simple,  fast to calculate, progressive,  fully deterministic, exhibits fast convergence rates (low variance), with desirable isotropic, blue noise spectra and therefore results in minimal artifacts.  So here goes… Figure 3. Rendering in 3D (left). Undesirable aliasing artifacts often arise in sampling processes (right). Source:[Perrier, 2018]   Outline simple jittered (finite) point sets jittered quasirandom point sets jittering of (infinite) quasirandom sequences QMC convergence comparisons summary example code Simple Jittered Sampling The most basic jittered (perturbed) point sets are the jittered grids. In this case, the sample space is subdivided into a grid of $m \times m$ square regions. Instead of merely selecting a point at the center of each cell, any (pseudo-) random point within each cell is selected.The advantages of jittered grid sampling is that it is trivial to understand,  simple to implement, and extremely fast to compute. It has the obvious limitation that $N$ must be a perfect square. Although rectangular stratified sampling can ease this limitation by simply requiring $N$ to be the product of two appropriate integers, I describe a different approach that is is informative, very elegant but less well known. Without loss of generality, assume that our sample space is the unit square $[0,1)^2$. We define the Fibonacci lattice as of $N$ points where: $$ \pmb{p}_i(x,y) = \left(  \frac{i-0.5}{N},  \left\{ \frac{i}{\phi}  \right\}\; \right); \quad \textrm{for } \; i=1,2,3,…, N. \tag{1}$$ where $\phi \simeq 1.618033988$ is the golden ratio, and $\{ x \} $ denotes the fractional part of $x$. Note that in computing, this fractional operator is usually written and implemented as $ x \; \%1$ (pronounced ‘mod 1’), and as this post as mainly geared towards computing applications we shall use this notation from here onwards. Note that the Fibonacci lattice is a very elegant method to place any number of points in the unit square in an extremely regular manner. For example, figure 2a shows the standard Fibonacci lattice for N=41 points. For jittered sampling, instead of selecting the center of each disk for our sample points, we can simply select a point $\tilde{p_i}$, inside each of the $N$ disks (see figure 4b). $$ \tilde{\pmb{p}_i}(x,y) =  \pmb{p}_i(x,y) + \pmb{\epsilon}_i (\delta_0/ \sqrt{N}),\quad \textrm{where } \delta \simeq 0.78 \tag{2}$$ where $\epsilon_i(r)$ denotes that for the $i$-th term, a random sample is drawn from a circle centered at zero with radius $r$. Note that a simple direct method to draw a sample uniformly from a circle of radius $r$ is to first draw two uniform random samples $u_1$ and $u_2$, each from the interval $[0,1)$ and then apply the following mapping: $$ \pmb{\epsilon}_i(r) = \left( r \sqrt{u_1}   \cos(2\pi u_2), r \sqrt{u_1}  \sin(2 \pi u_2) \right) \tag{3}$$ Of course, you can reduce the degree of jitter simply by selecting a smaller value of $\delta$.     The point distributions, approximate radial power spectra and the fourier spectra of these point sets is shown below.     If you do know $N$ in advance, then this is a simple but effective way to produce an isotropic point set. Interestingly, although intuitively there is a consistent distance between neighbouring points, the fourier spectrum of the jittered Fibonacci point set does not show the classic blue noise spectrum.  Nonetheless this jittered version is very good for many applications and contexts! However, there are known to be many excellent ways to produce isotropic point sets if $N$ is known in advance. This post mainly explores the lesser studied area of direct constructions for open/infinite isotropic blue noise point sequences, where $N$ is not known in advance. Furthermore, although the jittered grid and jittered Fibonacci lattice are extremely fast to compute, they also both suffer the disadvantage that the ordering of the points is ‘line-by-line’. Sampling methods where consecutive points are not proximal, but rather are more ‘evenly spread’ out across of the entire sample space, are often preferred. This property is often termed ‘progressive sampling’. Thus, we explore the viability of constructing progressive isotropic blue noise point sets and sequences. Jittering quasirandom point sets A low discrepancy sequence  quasirandom sequence is an (infinite) sequence points whose distribution is more evenly spaced than a simple random distribution but less regular than lattices. These sequences are useful in many applications, such as Monte Carlo integration. There are many types of quasirandom sequences, including Weyl/Kronecker, van der Corput / Halton, Niederreiter and Sobol.  I discussed these, including a new one in extensive detail in my previous post: “The unreasonable effectiveness of quasirandom sequences“. This post focuses on just one of these – the $R_2$ sequence. For comparative purposes, the Halton sequence and the Sobol’ sequence are also discussed in this post, as these are the most common in current QMC usage. A brief summary of each of the low discrepancy quasirandom sequences discussed in this post, is as follows. $\pmb{R_2}$: The $R_2$ sequence is a Kronecker/Weyl/Richtmyer recurrence sequence, introduced in my previous post. It is based on a generalization of the golden ratio, and so can be considered a natural way to construct a generalized Fibonacci lattice that wraps both left-to-right as well as top-to-bottom. It is defined as: $$ \pmb{t}_n = n \pmb{\alpha} \; \% 1,  \quad n=1,2,3,… \tag{4}$$ $$ \textrm{where} \quad \pmb{\alpha} =\left( \frac{1}{\varphi}, \frac{1}{\varphi^2} \right) \;\; \textrm{and} \; \varphi\ \textrm{is the unique positive root of }  x^3=x+1. $$ Note that $ \varphi = 1.324717957244746 \ldots$, often called the plastic constant, has many beautiful properties (see also here). Halton Sequence: the $d$-dimensional Halton sequence is simply $d$ independent Van der Coruput reverse-radix sequences concatenated in a component-wise manner. This post keeps to the common convention of selecting the first $d$ primes numbers for the bases. Thus, in this post, I always use the {2,3}-Halton sequence. Sobol Sequence: The Sobol’ sequence is known to have many excellent low discrepancy properties, including fast convergence rates. This post constructs the Sobol sequence (along with parameter selection) via Mathematica version 10. The Mathematica documentation states that the Sobol implementation is “closed, proprietary and fully optimized generators provided in Intel’s MKL library“. Unlike finite point sets, point sequences are infinite. A key advantage of infinite sequences, is that if the error after $N$ terms is still too large, one can simply sample additional points without having to recalculate all prior points. And an advantage of sampling based on quasirandom sequences is that they are progressive. That is, each successive point is proximally distant from the previous one, and the points fill the sample space evenly as the number of sample points increases. Figure 5 shows these three prototypical low discrepancy sequences, along with their fourier transforms and radial power spectra. We can see that: the $R_2$ points are very evenly spaced out. For many applications, this distribution is ideal, but for other applications it is ‘too even’ clear patterns are visible in all  the low discrepancy sequences, which often results in undesirable artifacts appearing in the QMC rendered images. the fourier spectra for all the low discrepancy sequences consists of set of spatially regular peaks none of the quasirandom point sets are isotropic the radial power spectrum for the $R_2$ sequence consists only of a discrete number of peaks. the radial power spectra for the Halton and Sobol sequences are continuous but not smooth. That is, rendering artifacts are a direct result of (i) discrete bright peaks in the fourier spectra, and/or (ii) discrete peaks in the radial power spectrum. A quick note on the term ‘blue noise’ may be useful here… In a loose sense, the $R_2$ distribution could be called blue noise because its low (blue) frequencies are clearly suppressed. However, it only consists of discrete peaks, and so may be better termed ‘discrete blue noise’ (or not called blue noise at all). In its more common context, ‘blue noise’ requires (i) an isotropic fourier spectra and (ii) a radial power spectrum that is flat like white noise, but with the low frequencies completely suppressed. For an excellent review of different blue noise spectra see [Heck, 2012]   Further, some authors [Perrier, 2018] note that most implementations of Poisson Disc sampling methods, such as Mitchell’s best candidate, are not strictly blue noise as their power spectra is small but not strictly zero at the lowest frequencies. In this sense, may be better termed ‘approximate blue noise’. They further comment that the difference between almost zero and exactly zero, can have a significantly detrimental effect on convergence rates.) We would like to combine the anti-aliasing advantages of jittered sampling with the fast convergence rates associated with low discrepancy quasirandom sequences, and thus: the main aim of this post is simple: to determine the minimum level of jitter to apply to the quasirandom points, such that the all fourier peaks are suppressed, the transform becomes isotropic, and the radial power spectrum for high frequencies becomes flat and even – ideally whilst keeping the lower frequencies suppressed. Intuitively this critical level of jitter would be closely related to the typical distance (divided by $1/\sqrt{N}$) between neighboring points. Figures 6a shows that the minimum separation between neighboring points for the $R_2$, Halton and Sobol sequences. Note that $R_2$ sequence is the only low discrepancy sequence where the minimum separation between neighboring points falls as slowly as $1/\sqrt{N}$. For Halton and Sobol the separation falls at much faster rate — closer to $1/N$. Similarly, figure 6b shows that the mean  separation in $R_2$ between neighboring points is at least $\simeq \delta_0/ \sqrt{N}$ for $\delta_0 /\simeq 0.76 $, all values of $N$.  These figures suggest that jittering $R_2$ within a radius of half this value would be ideal. Empirical tests verify this. Using this approach, the empirically determined critical radius for the jitter for $R_2$, Halton and Sobol equation are given in the following table. Their corresponding  point distributions are shown in figure 5. Note that the jitter for Sobol’ sequence has to fall much slower than $\sim … Continue reading A simple method to construct isotropic quasirandom
blue noise point sequences