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 rays coming from every direction that the object's surface reflects back to a given pixel. Unfortunately, a computer cannot compute an integrand directly. We therefore have two possibilities to solve this issue:

1. We find an analytical expression to remove the integrand (statistic based strategy).

2. We numerically approximate the equation by taking random samples in the integration domain and approximating the integrand value using Monte Carlo methods.

Here, we focused on numerical integration and sampling theory. Sampling is a fundamental part of numerical integration. A good sampler should generate points that cover the domain uniformly to prevent bias in the integration and, when used in Computer Graphics, the point set should not present any visible structure, otherwise this structure will appear as artefacts in the resulting image. Furthermore, 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 artefacts when used in Rendering.

The goal of this PhD was the study of the theoretical characteristics and practical impact of samplers combining both the Blue Noise and Low Discrepancy properties.

2. Approximer numériquement l'équation en tirant des échantillons aléatoires dans le domaine d'intégration et en en déduisant la valeur de l'intégrale via des méthodes dites de Monte Carlo.

La difficulté d'une représentation statistique de la surface réside dans le fait qu'une représentation statistique globale n'est pertinente que dans des cas très restreints (distribution gaussienne, pas de problèmes de visibilité, uniformité de l'objet ...).

Lors de cette thèse nous nous sommes donc plutôt intéressés à la deuxième solution, celle de l'intégration numérique, basée sur la théorie de l'échantillonnage. Ce sont ces méthodes qui constituent le noyau dur de cette thèse, et qui ont données lieu au développement de deux nouveaux échantillonneurs. Notre première méthode génère des échantillons 2-D ayant une excellente distribution spatiale. Elle est étendue dans notre deuxième méthode, qui peut produire des échantillons N-D ayant une distribution de qualité équivalente. Cette seconde méthode permet également de contrôler la quantité d'échantillons localement et en offre la possibilité de rajouter incrémentalement des échantillons tout en préservant la qualité de la distribution.

Dans la suite, nous présenterons tout d'abord la théorie de l'échantillonnage, avant de détailler nos contributions.

Théorie de l'échantillonnage

Dans le cadre du Rendu, si on ne peut pas calculer analytiquement l'intégrale, on va approximer sa valeur en simulant un certain nombre n de rayons lumineux. Cela revient donc à un problème d'échantillonnage; le problème étant de choisir correctement ces rayons pour avoir la meilleure approximation possible de l'intégrale. Pour ce faire, on va utiliser des méthodes dites de Monte Carlo et essayer d'approximer I grâce à l'équation (0.2).

E(I) = n i=1 p i f (x i ) (0.2)

Premier échantillonneur 2D

Ce premier échantillonneur était la première tentative pour concilier à la fois une discrépance basse et une distribution bruit bleu. En utilisant des ensembles de points 1-D connus pour être basse discrépance, on peut construire un ensemble 2-D où chaque ligne et chaque colonne est elle même un ensemble basse discrépance. Il a été prouvé qu'un tel ensemble 2-D est lui même un ensemble basse discrépance. Ensuite, en réarrangeant correctement cet ensemble 2-D, on peut changer la distribution des points. Cela nous permet d'obtenir une distribution de points bruit bleu (ou d'autres distributions), sans changer la discrépance de l'ensemble initial.

Pour réarranger l'ensemble de points, on permute les coordonnées x des points puis les coordonnées y. Intervertir des points formant un ensemble basse discrépance 1-D ne change pas la discrépance de l'ensemble, du coup notre ensemble 2-D résultant est toujours composé de lignes et de colonnes qui sont elles même basses discrépance. Ce qui fait que notre ensemble final est basse discrépance, tout en ayant une distribution de points différente. Pour obtenir une distribution ciblée, on précalcule et on stocke l'ensemble de permutations sur x et y à appliquer et l'échantillonneur n'a plus qu'à lire ces données et à faire la permutation. Tout ceci est illustré Cette méthode a le mérite d'avoir réussi à prouver que la basse discrépance et la distribution bruit bleu étaient parfaitement compatibles. Elle est cependant limitée dans ses applications pratiques. Elle ne génère que des points 2-D, n'est pas séquence, et ne permet pas de contrôler localement la densité des points.

Second échantillonneur N-D et Séquentiel

Ce second échantillonneur est lui aussi bruit bleu et basse discrépance. Il exploite le fait que si un ensemble est bruit bleu dans chaque sous partie (suffisamment grande) du domaine, il le sera également sur le domaine entier.

Notre méthode part donc d'une séquence basse discrépance qu'elle subdivise en tuiles carrées. Ensuite, on applique nos permutations locales, sur chaque tuile. Nous avons prouvé formellement dans le papier que ces permutations nous permettent de préserver la discrépance de l'ensemble entier malgré le fait qu'elles soient locales. Parmi toutes les permutations locales possibles, nous avons précalculé celles qui correspondaient au mieux à une distribution bruit bleu, et nous les avons stockées.

De plus, grâce à des permutations globales, dont nous avons également prouvé qu'elles préservent la discrépance, nous pouvons rendre notre échantillonneur séquence. Cela signifie que si nous avons un ensemble de points répartis dans des tuiles à un niveau l, nous pouvons resubdiviser ces tuiles au niveau l + 1 pour rajouter des points tout en préservant la position des points déjà présents.

Enfin, comme notre système est basé sur des tuiles, on peut choisir de ne subdiviser qu'une partie des tuiles pour changer localement la densité des points. empty page viii

Acknowledgement

Tout d'abord, je m'excuse auprès de tout ceux qui ont (ou qui vont) lire ce document ... Un grand merci aux rapporteurs et membres du jury qui eux n'ont pas eu le choix ! Pour le reste je suis pas sûre de savoir quoi marquer ici, j'ai toujours été nulle avec les mots ... Je n'ai rien contre faire vos éloges individuels, mais vous êtes tous tellement géniaux que ça augmenterai massivement le nombre de pages de ce document et après j'aurai toutes les associations contre la déforestation sur le dos. Je ne tiens pas non plus à faire une liste de noms parce qu'on fini toujours par oublier quelqu'un et c'est super vexant. Du coup je veux juste dire à tous les potos de labos (oui David et Victor vous êtes aussi mes potos), à tous mes amis et à toute ma famille à quel point je vous aime fort. Sans vous ces 3 années n'auraient été qu'un long cauchemard, et j'aurai probablement tout laissé tomber au bout d'une seule année sans regrets. Grâce à vous non seulement j'ai tenu suffisamment longtemps pour écrire ces lignes, mais maintenant que c'est fini, je pars en plus avec un gros pincement au coeur en pensant que je ne vous verrai plus (ou plus assez à mon gout !). Un grand merci pour tout ce que vous m'avez apporté (burger, body combat, équitation, jeux de société, séances de bitch/gateaux, robot de cuisine trop bien, et j'en passe ...) Je finirai avec cette citation tirée d'un chef d'oeuvre cinematographique moderne

Friends have a way of making even the worst of times into something pretty great. Nowadays, computer generated images are ubiquitous. They appear in movies, commercials, video games and have become almost undistinguishable from real pictures. This "realisticness" is partly due to the modern technologies. They increase so drastically our processing power, that it is now possible to try simulating the physical light field in the scene. The issue is that a physical light field in the real world is something that is infinite and continuous. Computers cannot handle any of those properties, and therefore, to simulate this light field, we need to rely strongly on sampling and discretization.

To generate a synthetic image we need to first be able to sample this simulated light field (meaning being able to measure the amount of light in a particular point of the scene), and second, we need to be able to reconstruct a 2-D image from those sampled values. Our eyes do this reconstruction thanks to a set of cells covering our retinas. Each of those cells receives a discrete portion of the real light field in front of us. Those values are then sent to the brain that reconstructs an image from it. To create a synthetic image from a simulated light field, we replicate this process of sampling a set of discrete values from the simulated light field and use tools such as convolutions or Fourier transform to reconstruct an image from it. Interested readers can find a detailed description of this process in appendix A on page 168.

Our main focus will then be on how to simulate and sample this light field. As every object in a scene receives light and reflects light, evaluating the amount of light that reaches a particular point in a scene is very computationally expensive. It requires evaluating all the light incoming to this point from every other point on every possible objects in the scene. Mathematically, this translates as a recursive integration, which cannot be solved analytically. It requires to be approximated, using numerical integrations methods, which strongly relies on sampling. In this document, we will focus on the use of Monte Carlo estimators as our tool to solve this numerical integration. Such light field simulation will be described Section 1.1, with a focus on Monte Carlo estimators in Section 1.1.2 on page 4.

The Rendering process

To render an image, we need to simulate the flow of light bouncing within a scene. This implies recursively integrating for every single point of the scene the light incoming from every possible direction. To solve this set of complex integrands, we use Monte Carlo estimators, that rely on sampling to approximate the result of an integration. In the following, we first detail how we simulate the flow of light and the different integrands that need to be accounted for. Then, we see how we can use Monte Carlo estimators to approximate those integrands. i m a g e p l a n e

x p

The light field in the scene is represented as a set of discrete light rays. The amount of light within each ray is computed by sending this ray from the image plane towards the scene, and bouncing it against the geometry until it reaches the light source.

Path-Based rendering

Generating a synthetic image first requires a description of the scene to be rendered. This description contains a set of objects present in the scene (often referred to as the geometry of the scene), each of them being associated with a material, a color, a texture, a motion, etc. A scene description also contains a set of light sources and a description of the virtual photo camera taking the picture (its position in the scene, its lens and focal, its looking-at point etc.). This camera defines a virtual image plane in the scene, usually discretized into pixels.

To render a scene, we compute the amount of light incoming from the scene to each pixel of the image plane. This set of discrete radiance values will then be the input for a reconstruction process which will generate the final image (as described in Appendix A on page 168). This radiance can be seen as the amount of light reaching a point x from every possible directions.

To evaluate this, we use path tracing algorithms. As this is a very active field of study and as it is a bit out of scope of this document, we will only present the simple naive algorithm here (Fig. 1.1). Interested readers can refer to [T+12; Pha+16; Vea97] for more details.

Path tracing algorithms [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF] aim at solving the following equation, first devised by Kajiya in 1986, I(x) = L e (x) + S ρ(p, x)I(p)dp.

(1.1)

Given here in its simplest form, this equation, known as the Rendering Equation, evaluates the amount of light I that reaches a point x on the image plane. To do so, it sums the intensity of the emitted light of x with the amount of light ρ(p, x) that each point p of the whole geometry surface S reflects towards to point x. One of the difficulties with this equation is that ρ(p, x)

is not a continuous function, it includes a visibility term, which is a binary factor with values being either 1 or 0, that accounts for the visibility discontinuities within the scene.

Note that light is composed as a set of wavelengths, which determines its color. Since the colors, materials, or geometry of the objects in the scene can affect each wavelength differently, this equation needs to be solved for each of them. However, this problem is out of the scope of this document and, to simplify the notation, we will assume that we are always working with a single given wavelength.

The issue that interests us most here, is that Equation ( 1.1 on the preceding page) relies on recursive evaluation of complex integrands. As we assume the scene to be a continuous signal, those integrands translate as infinite sums. It is thus impossible to evaluate them directly (or only in very specific cases [START_REF] Bruneton | A Survey of Non-linear Pre-filtering Methods for Efficient and Accurate Surface Shading[END_REF]). Instead, we will rely on sampling-based numerical integration methods. More specifically, we will use Monte Carlo estimators which will be described in details in section 1.1.2 on the following page. For now, we just describe Monte Carlo estimators as a way to approximate the integrand of a function by evaluating N points of this function, and averaging them. Following this, Equation 1.1 on the preceding page becomes

I(x) ≈ L e (x) + N i=1 ρ(p i , x)I(p i ) (1.2)
where p i is the i th sample. The quality of the approximation is now dependent on the error from the Monte Carlo estimators, which is related to the way the samples p i are chosen (see Section 1.1.2 on the following page). Note that in this whole document, we only focus on unbiased estimators, meaning estimators that will converge to the true value of the integrand when using a infinite number of samples. Therefore we will often implicitly ignore bias. Note that other estimators could have been used, but we focus on estimators that do not assume prior knowledge of the function to integrate. For more details please refer to [Gla95;[START_REF] Hammersley | Monte Carlo Methods[END_REF].

To evaluate each I(p i ), we combine Monte Carlo estimators with Path tracing algorithms. This means that we represent our light field as a set of discrete light rays. Then, we evaluate the light field value for each ray by simulating how it interacts with the scene.In the real world, rays are coming from the light before reaching our eyes. In this algorithm, we do the opposite. The rays are thrown from the image plane and bounce within the scene either a finite amount of times or until they reach the light.

For each pixel of the image plane, we will sample a certain amount of points {p i }. Each p i will be used as the starting point for a light ray. Each ray is then tested against the scene to compute the intersection point p i between this ray and the geometry. Once this intersection point is found, we sample a bouncing direction of the ray and we compute the amount of light that bounces in this direction. To do so we use a Bidirectional Reflectance Distribution Function (BRDF). This function takes as input an incoming light direction (the direction of the ray) and an outgoing light direction (the sampled direction) and returns the amount of light that was reflected. This BRDF function is specific to each possible material of each object of the scene. Then, we once again test the bounced ray against the scene to find its new intersection point with the geometry, sample a bouncing direction, compute the reflected light, etc. This recursion ends when the ray hits a light source, or when a maximal number of recursion is achieved (Fig. 1.1 on page 2).

Following this recursive algorithm, the evaluation of a single sample can become quite costly. The challenge is thus to have an approximation as good a possible with as little samples as possible. This is the focus of the next section, along with the theory of Monte Carlo estimators.

Monte Carlo estimators

Monte Carlo estimators are used to numerically approximate a complex integrand. They rely on a set of samples, that are used to evaluate the function to integrate. Then, the values of this function for each sample are averaged to give an approximation of the integration.

More formally, to use classic Monte Carlo estimators, we first need to generate n independent, uniformly distributed random variables {x (i) } n in Ω, that we will refer to as samples. Then, an estimation of the mean value of a function f can be obtained via the unbiased estimator:

I n := 1 n n i=1 f (x (i) ). (1.3)
Note that this function f is defined on a s-dimensional domain Ω. Thus, all the samples x (i) must be generated within this same domain (following [Gla95] for notations). In most cases, Ω = [0, 1) s but this is not mandatory.

The integration error Δ when using such an estimator is defined as

Δ := |I n -I| (1.4)
where I := Ω f (x)dx is the true integration value. Note that this error depends on the sampling pattern used, the number of samples used, and the integrand function.

For stochastic samplers, we often use instead of this error the Mean Square Error (MSE), that will be defined just a bit below. First, we denote

Bias := |I -I n | (1.5)
and V ariance := (I n ) 2 -I n 2

(1. 6) where . represents the expectation of a random variable.

From [START_REF] Durand | A frequency analysis of Monte-Carlo and other numerical integration schemes[END_REF][START_REF] Subr | Fourier analysis of stochastic sampling strategies for assessing bias and variance in integration[END_REF], the MSE can then be expressed as MSE(I n ) = Bias 2 + Variance.

(1. 7) In this thesis we only work with unbiased stochastic point sets, our MSE depends thus entirely on the variance of I n . For classical Monte Carlo estimator with pure random samples, this variance is

σ 2 := σ f n , (1.8)
where σ f is the variance of the integrated function f . As the variance is function of the number of samples, the more samples we use the less error we observe (for unbiased point sets).

With a MSE in O(1/n), the error in integration Δ is in O(1/ √ n). Thus, to divide this error by 2, we need to use 4 times more samples, which is rather unsatisfying. Furthermore, even though when taking a infinite number of samples in a pure random way we are guaranteed that we will cover the whole domain uniformly, there is no guarantee, if unlucky, that when taking a given number n of samples, they will not all clump into a small portion of space. For all those reasons, using pure random samples is close to be the worst scenario for unbiased Monte Carlo integration. Therefore, several improvements were developed, reducing the variance of this estimator by alleviating this clumping.

One of the most popular is the use of stratified sampling (Section 3.1 on page 36). This sampling subdivides the domain into k equal regions and generates random samples within each region. This sampling has provably a better variance than pure random sampling [START_REF] Hammersley | Monte Carlo Methods[END_REF], proportional to the number of regions k and the way they are chosen. The exact variance function when using stratified sampling is quite complex and thus we will not give it here. Interested readers can find it in [Gla95], Section 7.4.3, Equation (7.35).

Another improvement, which interests us more here, is the usage of samples that are no longer randomly chosen. By choosing such samples in order to guarantee an optimal covering of the domain, you may reduce the error of the approximation. This category of estimators are called Quasi Monte Carlo (QMC). We point that such estimators were not originally designed for computer graphics applications but for financial and simulations contexts.

When using deterministic sets of samples, we can no longer rely on the MSE value, as the approximation I n is not longer a random variable. The quality of QMC estimators is thus only measured from their error, which is

Δ = k log(n) b n (1.9)
for constants b and k [Gla95].

However, even though such samplers lead to a lower error and a better approximation for the color of a single pixel, this is not enough to ensure a better final image. An image is composed from a whole set of pixels, and therefore from a whole set of integrands, all solved independently. The way those integrands are correlated or not also impacts the final rendering.

Sampling for Monte Carlo Rendering

There are two specificities when using Monte Carlo estimators for rendering. The first one is that neighbouring pixels should be decorrelated, to avoid structural artefacts in the final image (Fig. 1.2 on the next page). The second is that the integrand to solve is recursive, therefore, even for a single pixel, we need to solve several integrands all related to each other. It is interesting to note that in most situation, we are working with recursive 2-D integrands, since we are mostly sampling directions or points on a surface.

The structural artefacts appear mostly when using QMC estimators for rendering. QMC sampling patterns were originally devised for applications in finance and simulation. They have been proven to lead to a very low integration error, along with being easy to implement and computationally efficient, even in very high dimensions. Those properties are very often due to the analytical definition of such samplers. However, such analytical definition often implies that they present a recursive and regular spatial organisation of samples. This structure correlates the integration error over neighbouring pixel, which leads to disturbing structures in the final image. This is studied in more details in Section 2.1 on page 14 and in Section 3.3 on page 57. Instead, when using random samplers, this structure turns into noise [START_REF] Mark | Antialiasing Through Stochastic Sampling[END_REF], which is directly related to the variance in the estimator. It can never be completely removed but it can become imperceivable to the human eye, as the variance reduces (Fig. 1.2 on the facing page). Roughly, using QMC estimators leads to higher artefacts but a faster convergence and MC estimators presents a slower convergence but with lesser artefacts.

The second specificity of Monte Carlo Rendering is that even though we are solving a recursive multi dimensional integrand, not all dimensions are directly correlated. Assume that we are sampling a light ray that has N bounces. We will need one 2-D sample for the light plane. Then, we will need N 2-D sample for the bouncing directions. We can also need a 1-D sample for the motion blue, or another 2-D sample for the lens. The point is that we do not necessarily need a s-D sampling pattern. One can also use several 2-D sampling patterns, while being careful to decorrelate them and remove any structure. However, when using N 2-D sampling patterns, this decorrelation is simulated through random shuffling of samples, which harms the quality of the approximation of the integrand. This will be detailed in more specifics in Chapter 4 on page 74.

Contributions

In this PhD, we focused on developing proofs and algorithms to generate sets of s-D samples (also referred to as point sets, distributions of samples, or sampling patterns) that will be as optimal as possible regarding Monte Carlo rendering.

We tried to address 3 major issues, • Removing the bias in integration by guaranteeing that our point sets are uniformly distributed (Figure 1.2 on the preceding page (c)).

• Removing any structure in our spatial distributions to remove all aliasing in the final image (Figure 1.2 on the previous page (d)).

• Finally, reducing as much as possible the variance in Monte Carlo/Quasi Monte Carlo estimation (Figure 1.2 on the preceding page, (b)).

To assert all those issues, we tried to combine the uniformity and efficiency properties of QMC samplers with the irregular samples placement of stochastic point sets. Joining those properties apply additional constraints on the spatial samples organisation, which were though to be unreconcilable. We showed that they are in fact reconcilable, by providing several algorithms generating such point sets, along with their formal proofs. Some possible naive approaches to generate such samplers are presented in chapter 4 on page 74 and our main contributions are presented in chapter 5 on page 93 and in chapter 6 on page 107.

As a preliminary, Section 1.4 will present a few basic notations and definitions regarding point processes. Those definitions will be further developed in Chapter 2 on page 14 where we will see in deeper details how the quality of a sampling pattern regarding Monte Carlo Rendering can be assessed.

Basic definitions

When using a MC or QMC estimator, we need to generate a set of variables x (i) . We denote point set or point pattern a set P n of n indexed s-D points P n = {x (1) , x (2) ,

• • • , x (n) }
where all x (i) ∈ Ω, with Ω ⊂ R s (Fig. 1.3 on the facing page). Note that as P n is a set, we implicitly assume that there are no duplicates in this set, more formally, we assume that ∀i, j, x (i) = x (j) . Point sets are generated by samplers, which we will denote S.

We denote S(Ω) all point sets generated from the sampler S in the domain Ω. Note however, that in most rendering situations, the samples are defined in the unit hypercube [0, 1) s , before being rescaled to the integration domain (which may be the surface of the light, a pixel, a direction hemisphere, ...). This is why all along this document, we will ease the notations by implicitly assuming that Ω is the unit hypercube (formally defined as Ω = [0, 1) s ), unless specified otherwise.

Deterministic samplers generate from a given set of parameters a unique set P n (note that the value n is a parameter of S). Typically, such samplers are used with QMC estimators. On the contrary, stochastic samplers can generate several point sets from a given set of parameters. A point set P n generated from a stochastic sampler is referred to as a realisation of this sampler.

We assume that if we denote A the set of point patterns that present a given characteristics, we have S(Ω) ∈ A, meaning that a stochastic sampler generates different sets but always with similar characteristics (see Chapter 2 on page 14 for a list of criteria that may characterize a sampling pattern). Such sets are typically used for MC estimators.

In her book, Illian [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF] presents a thorough statistical analysis for samplers (which she refers to as point processes). We note that when S is stochastic, S(Ω) is a random variable. For simplicity, we will assume that this is also true for deterministic samplers, with all realisations of S(Ω) being equal for a given set of parameters. We can then differentiate different types of point processes. A point process can be finite, infinite, isotropic, stationary and/or ergodic.

Finite/Infinite point processes A sampler S that can generate any n number of samples over any domain Ω is an infinite point process. Finite point processes are point processes defined over a single given Ω. They are mostly related to real world situations, for example, the repartitions of the seeds around a tree, or the number of pores in a Swiss cheese.

P(S(Ω)) = P(S(Ωx))

(1.10)

and isotropic if ∀α ∈ [0, 2π] R α S(Ω) = {R α x (1) , R α x (2) , • • • , R α x (n) } (1.

Basic definitions

Density of a point processes We also define the density of a point set as the average number of samples inside a region R of volume V R around a sample x.

λ(x) = R(x) V R (1.13)
This density is constant for isotropic and stationary point processes. A sampler generating sets with a non constant density is sometimes called a non-uniform sampler.

Ergodic point processes

The ergodicity of a sampler S refers to samplers that are so statistically invariant that the analysis of a single realisation is statistically meaningful. We point out that a deterministic sampler is necessarily ergodic.

Sequentiality We will end this section by defining the notion of sequentiality for a sampler. A sampler is said to be a sequential sampler if it does not require a prior knowledge of the number of samples to generate. Such samplers are able to iteratively add n new samples within an already existing set, while preserving the characteristics of the set. We denote sequence a sampler that is both sequential and present highly uniformly distributed samples (low discrepancy). We will see how we can characterize this uniformity in Chapter 2 on page 14.

Finally, as many of the samplers we will study rely on binary manipulations over positive integers, we define as

x = N k=0 a k b k (1.14)
the representation of a integer x ∈ N in base b, where a k ∈ Z/Z. Each a k value will be denoted a digit of x, or a bit of x if b = 2. The value N is the number of digits needed to encode x in base b. Note that each a k value belongs to the domain Z/bZ. For simplicity, we will also denote x from the list of digits of its b-ary representation as

x = (a N , • • • , a 0 ) b . (1.15)
Similarly, if a value x belongs to the set [0, 1), it will be denoted as

x = N k=0 a k b -k-1 , (1.16) and x = (.a 0 , • • • , a N ) b .
(1.17)

Finally, we will denote (a N , • • • , a 0 ) or (a 0 , • • • , a N ) the vectors made from the list of digits of a value x in base b.

Note that any value x ∈ N can be turned into its equivalent y ∈ [0, 1) (and vice versa) with the following operations

y = x/b N (1.18) x = y • b N . (1.19)
This means that all those notations are roughly equivalent, and in this document, when there is no ambiguity, we will often implicitly switch between them. Table 1.1 on the next page summarizes all the notations that will be used in this document.

In the rest of this document, we will rely on those definitions to present the different characteristics of a set of samples that impact the quality of the final rendered image (Chapter 2 on page 14). We will then present existing samplers (Chapter The total number of dimensions x A point in s-D x (i) The i th sample of a set P n x i

The i th coordinate of x x (i) j

The j th coordinate of

x (i) P n A point set of size n Ω
The domain of definition of S S(Ω)

A point set generated from the sampler S Chapter 2

F f
The Fourier transform of a function f

F -1 f The inverse Fourier transform of a spectrum f P f The Power spectrum of f D * The l ∞ discrepancy D The l ∞ extreme discrepancy D o The l ∞ isotropic discrepancy T 2
The l 2 discrepancy D 2

The generalized l 2 discrepancy Chapter 3

(a N • • • a 0 ) b The b-ary representation on a value x ∈ N (.a 0 • • • a N ) b The b-ary representation on a value x ∈ [0, 1) (a N • • • a 0 )
A vector of values with a k ∈ Z/bZ

Characterizing a set of samples

There are several ways to characterize a point set. Depending on the application, one can focus on visual quality, stochasticity, ease to control, etc. In our case, we will use those sets for Monte Carlo rendering. We will thus characterize our samplers with characteristics related to the integration error of their MC estimator. This leads us to three main characteristics. First, the uniformity of the point set (Section 2.1), second, its Fourier spectrum, (Section 2.2 on page 20) and third, its spatial organisation (Section 2.3 on page 25). We chose those measures specifically, since they are known to be related to the integration error and integration variance (Section 2.4 on page 28), and therefore are related to the expected quality of the final rendered image.

Uniformity

A set of samples can be characterized by how close it is from the perfect uniform distribution.

In other words, how well the samples are spread out over a given domain. There are several ways to measure this uniformity.

Some tests rely on a supposed lattice structure of the distribution, like the spectral test [START_REF] Coveyou | Fourier analysis of uniform random number generators[END_REF]. This test, defined only for lattices, measures the largest distance d s between adjacent parallel hyperplanes that together cover the whole point set. The smaller this distance, the better the uniformity of the lattice (Fig. 2.1). Other tests apply on point sets quantized on a grid of a given resolution (Figure 2.2 on the next page left), namely the test of Congruent Boxes (Figure 2.2 on the facing page right). Roughly, the domain is partitioned according to a grid in base b, and if each box of this grid contains the same number of samples, the point set is considered uniform [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF]. This test is used to determine if a point set is a digital net. Such nets will be of importance and will be presented in more details in Section 3.3.2 on page 59. The most common measurement of uniformity is the discrepancy of a point set. One of the main asset of this measure is that it is not specific to a type of point set. It quantifies the distance between a distribution of samples and the uniform distribution. In 1-D, we compute the distances between the Cumulative Distribution Functions (CDF) of the uniform distribution, and the CDF of a given 1-D point set P n , both defined in the domain [a, b). Note that in most situations, we have a = 0 and b = 1. The CDF of a point set (Fig. 2.3 on the following page) is defined as

d 2 =0.196 d 2 =0.128

Quantized Point set Subset of Congruent Boxes

F n (u) := 1 n Card({x, x ∈ P n , x < u}) (2.1) with u ∈ [a, b).
The uniform distribution, denoted U is defined as the distribution where each possible random variable is equally probable. A 1-D illustration of the distribution U(a, b), which is the uniform distribution of all the values in [a, b), is given Figure 2.4 on the next page.

We quantify the l ∞ distance between the distribution of a given point set to U using the Kolmogorov-Smirnov test. The distance N (P n ) between a point set P n and the uniform distribution is then

N (P n ) := sup u |F n (u) -F (u)| (2.2)
where sup u returns the maximum distance and where F (u) is the CDF of the uniform distribution (Fig. 2.5 on the following page). Note that we could also quantify similarly the l p distance D p (P n ) between the distribution of our point set and U,

N (P n ) := x |F n (x) -F (x)| p . (2.3)
Those definitions apply to 1-D point sets. We will see in the following sections the definition of the l ∞ star discrepancy (also referred to as star discrepancy [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]) for 2-D point sets, and the l 2 star discrepancy for s-D point sets.

l ∞ discrepancy

One of the most common of all the existing discrepancy measure is the l ∞ star discrepancy. This measure considers all axis-parallel hyper-rectangles B(v) in the domain [0, 1) s with a corner at the origin of the domain and the other one at position v (Figure 2.6, left). We call such hyper rectangles anchored boxes. We then define the function

α(P n , v) := Card(P n ∩ B(v)) n (2.4)
that counts the number of samples that fall inside an hyper-rectangle and expresses it as a fraction of the total number of samples n. We then compare it to the true area of

B(v), v 1 • • • v s .
This leads us to the final formula for the l ∞ star discrepancy,

D * (P n ) := sup v∈[0,1) s |v 1 • • • v s -α(P n , v)|. (2.5)
One can note the existence of other discrepancy measures than the l ∞ star discrepancy. First, the l ∞ extreme discrepancy that computes the discrepancy using unanchored axis-parallel hyper rectangles (Figure 2.6, middle). Note that the following relationship holds between the l ∞ star discrepancy D * (P n ) and the l ∞ extreme discrepancy D(P n ) [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF]:

D * (P n ) ≤ D(P n ) ≤ 2 s D * (P n ). (2.6)
Finally, it is worth noting the existence of the l ∞ isotropic discrepancy that computes the discrepancy using any convex shape (Figure 2.6, right). The following relationship holds between the l ∞ isotropic discrepancy D o (P n ) and the l ∞ extreme discrepancy [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF]:

D(P n ) ≤ D o (P n ) ≤ (4s √ s + 1)D(P n ) 1 s . (2.7)
As the extreme and isotropic discrepancies are extremely difficult to compute and are both related to the star discrepancy, we will not discuss them in details here.

Measuring the discrepancy of a point set characterizes the uniformity of this single set. However, a careful reader will have noted that, aside from pathological samplers, the discrepancy of a set decreases as the size of the point set increases. We can use this property to characterizes a sampler with the "speed" at which its discrepancy decreases. For example, for a pure random sampler in 2-D, the l ∞ star discrepancy decreases at a speed of 1 √ n , when for the regular grid, it decreases in 1

-(1 -1 n ) 2 [Lem09; Nie92].
The following lower bound on the l ∞ star discrepancy has been shown for s = 2 [START_REF] Schmidt | Irregularities of distribution[END_REF]:

D * (P n ) ≥ C log(n) s-1 n , (2.8)
where C is a constant dependant on s but independent of n [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF].

If you have a sampler with a discrepancy decreasing in O log(n) s-1 n as the number of samples n increases, this sampler is said to be a low discrepancy sampler. If this sampler is also sequential, we refer to it as a sequence. The l ∞ star discrepancy is defined using axis-parallel hyper rectangles anchored at the origin (Figure 2.6 on the preceding page, left). This means that when two samples have the same coordinate on either axis, a box of infinitely small size can contain several samples. This almost guarantees to break the low discrepancy property. Therefore, one of the first conditions (necessary but not sufficient) for a point set to be low discrepancy is that it never presents several samples aligned on either dimensions.

Computing the l ∞ star discrepancy can be done using explicit formulas for s = 1

D * (P n ) = 1 2n + max 1≤i≤n x (i) - 2i -1 2n (2.9)
and for s = 2, as detailed in [START_REF] Niederreiter | Discrepancy and convex programming[END_REF] (note that similar formulas exist for the extreme discrepancy).

It has been observed that the set of hyper-rectangles to consider can be discretized (Fig. 2.7 on the next page). Following the notations in [START_REF] Gnewuch | A randomized algorithm based on threshold accepting to approximate the star discrepancy[END_REF], we denote

X j (P n ) := {x (i) j , i ∈ {1, • • • , n}},
the set containing all the j th coordinates of a set P n . We also denote

X(P n ) := X 1 × • • • × X s ,
the set of all the points that can be made from all the coordinates of all the samples of P n . Now we can redefine the l ∞ star discrepancy as

D * (P n ) = sup v∈X(Pn) |v 1 • • • v s -α(P n , v)| (2.10)
which can be computed exactly using an enumeration algorithm in O(n s ) [START_REF] Niederreiter | Discrepancy and convex programming[END_REF] (Fig. 2.7 on the facing page). We stress that such an algorithm is only practicable for small values of n and/or s. The computation of the star discrepancy has been proven to be a NP-hard problem

v Γ 0 (P n ) Γ 1 (P n ) Fig. 2.7.: Computing the discrete l ∞ star discrepancy in 2-D.
in higher dimensions [START_REF] Gnewuch | Finding optimal volume subintervals with k points and computing the star discrepancy are NPhard[END_REF]. Note that the best known algorithm to compute the l ∞ star discrepancy is in O(n 1+s/2 ) [START_REF] Dobkin | Computing the Discrepancy with Applications to Supersampling Patterns[END_REF] and that several algorithms exist to approximate it in

s-D [Thi01; Gne+11; Sha10].
This discrepancy is not the only defined discrepancy measure, one can also compute l p star discrepancies. We find that the l 2 star discrepancy can be computed exactly in O(sn 2 ), which is useful as s increases.

l 2 discrepancy

The l 2 star discrepancy, denoted T 2 , is a variation of the l ∞ star discrepancy. Instead of using the l ∞ norm (sup), we use the l 2 norm :

T 2 (P n ) := [0,1) s (|v 1 • • • v s -α(P n , v)|) 2 dv 1 2 . (2.11)
It should be used carefully as it is known to both put an emphasis on points close to 0 and to be irrelevant for sets that contains too few samples relatively to the number of dimensions [START_REF] Matoušek | On the L2-discrepancy for Anchored Boxes[END_REF]. However, it can be computed very easily following Warnock's formula [START_REF] Turner | Computational Investigations of Low-discrepancy Point-sets[END_REF],

(T 2 (P n )) 2 = 1 n 2 n i=1 n j=1 s k=1 1 -max x (i) k , x (j) k - 2 -s+1 n n i=1 s k=1 1 -x (i) k 2 + 3 -s , (2.12)
where

x (i)
k is the k th coordinate of the i th sample in the set. This gives us an algorithm in O(sn 2 ) to compute T 2 .

Several versions on the l 2 star discrepancy were defined [START_REF] Hickernell | A generalized discrepancy and quadrature error bound[END_REF], each of them with a formula similar to Warnock's. Among them, the one that interests us here is the generalized l 2 star discrepancy, defined by Hickernell [START_REF] Hickernell | A generalized discrepancy and quadrature error bound[END_REF] with the following formula

(D 2 (P n )) 2 = 4 3 s - 2 n n i=1 s j=1 3 -x (i) j 2 2 + 1 n 2 n i=1 n i =1 s j=1 2 -max x (i) j , x (i ) j (2.13)
This value is of special interest since it is related to the integration error, as will be discussed in section 2.4.1 on page 28. The generalized l 2 star discrepancy can also be computed in O(sn 2 ). This makes it much more efficient in higher dimensions than the l ∞ star discrepancy, that has a complexity in O(n 1+s/2 ) [START_REF] Dobkin | Computing the Discrepancy with Applications to Supersampling Patterns[END_REF].

For all those reasons, when evaluating the uniformity of higher dimensional pointsets (in practice, for all sets with s > 2), we rely on the generalized l 2 star discrepancy.

Those discrepancies measures are of uttermost importance in Monte Carlo integration. If a sampler is not uniform, meaning if it does not verify

lim n→∞ D * (P n ) = 0,
it is not guaranteed to converge to the correct integrand value when used with Monte Carlo estimators. However, even if it is uniformly distributed, having a structured sampling pattern may lead to correlation between neighbouring samples, which leads to aliasing in the final renderings. To characterize the amount of aliasing generated from using a particular sampling pattern in rendering, we rely on spectral and spatial analysis (Section 2.2, and Section 2.3 on page 25).

Spectral analysis

A sampler can also be characterized by how the samples are spatially organized, and more specifically, by which frequencies arise in this distribution. First we will see how to compute the frequency representation of a signal and then, we will see how we can characterize a point set from this representation. Note that in [START_REF] Newbold | The Fourier transform and its applications[END_REF], interested readers can find a much more detailed description of all the tools presented here.

The Fourier transform

There are two mains representations for a signal. The most common is the spatial representation, which associate to a value x ∈ Ω a value y ∈ f (x), with x being the time or a position in space, and y the function value. However, it appears that many real-valued signals can be represented by a combination of sinusoids with different phases, amplitudes and frequencies. Such representation of a signal is called its frequency, or spectral representation. The transfor-mation from the spatial representation to the frequency representation of a signal is given by the Fourier Transform.

The representation of signals as a set of sinusoids is referred to as the Fourier series, defined in 1-D by

f (x) = a 0 + ∞ n=1 a n cos(xnω) + ∞ n=1 b n sin(xnω) (2.14)
where the a n and b n , that control the amplitude of the sinusoids, are defined in the complex domain C and where nω is the frequency of the sinusoids. Thus, the Fourier series can also be expressed as a sum of complex exponentials

f (x) = ∞ n=0
c n e ixnω .

(2.15)

Those Fourier series assume that the represented signal is periodic. For non periodic signals, we will rather use the Fourier transform. This operation translates from the spatial to the frequency representation of f . This transform associates to a frequency ω its amplitude and phase in the Fourier serie decomposition of f . It is given by the following equation in 1-D,

F f (ω) = R f (x)e -ixω dx (2.16)
for a function f (x), defined over R, with f presenting frequencies over the domain Ω. The real part of F f (ω) represents the amplitude of the sinusoid and its imaginary part the phase of the signal f at frequency ω. Note that

F f (0) = R f (x)dx.
(2.17) Figure 2.8 on the following page illustrates the relationship between the Fourier transform and the Fourier serie of a signal. In the following, we will often refer to F f (ω) as the spectrum of a signal f .

Note that the Inverse Fourier transform can transform back the spectrum of f into the original signal f as

f (x) = C F f (ω)e ixω dω (2.18)
The existence of this function allows to say that the frequency and spatial representations of a signal are dual.

For simplicity, we only gave formulas and examples for the 1-D Fourier transform. However, the Fourier transform is defined in s dimensions. A s-D signal presents a s-D spectrum, computed with similar equations. You just replace the scalar products between x and ω with s-D dot products in Equations ( 2.16) and 2.18.

However, we point out that for practical applications, we handle signals that are defined over a closed domain Ω. Therefore, to apply a Fourier transform to our point sets, we implicitly assume that they repeat periodically (to avoid oscillation in the Fourier transform due to border effects). Furthermore, we handle only a discrete representation of our signals, therefore, in practice, we perform discrete and periodic Fourier transforms, which are equivalent to Fourier series.

F 1 = 4 π sin(x) F F 0 -1 -0.
F 3 = 4 π sin(x) + 4 3π sin(3x) + 4 5π sin(5x) F F 3
F N = 4 π N n=1 1 n sin(n × x) F F N

The Power Spectrum

Here, we want to characterize our point sets by the amount of energy of each frequencies.

Therefore, what interests us is the amplitude of each sinusoids in the decomposition of our s-D point set P n . We denote by P Pn the Power spectrum of a point set. P Pn is formally defined as

P Pn (ω) = |F Pn (ω)| 2 . (2.19)
Using the norm of F Pn (ω) removes the imaginary part and preserves only the squared amplitude. This results in having all the values of P Pn in R, with the removal of the phase information that was contained in the imaginary part of F Pn .

A s-D point set, leads to a s-D power spectrum. So, for simplicity, we will often use its radial representation instead. This representation averages all the values of the Power spectrum between two s-D spheres, turning a s-D representation to a 1-D representation of our Power spectrum. This radial average is only truly valid for isotropic spectra. Still, for anisotropic spectra, it can still be used to identify peaks in the Power spectrum, a.k.a. frequencies that are so present that they are almost unaffected by averaging. All of this is illustrated Figure 2.9 on the next page.

Differential Analysis

Differential analysis is a tool was presented in [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF], and that can be thought of as a reformulation of the classical spectral analysis. It lets us compute the spectral properties of sets of uniform and non-uniform density by reformulating the spectral analysis as a difference between pairs of samples.

The Fourier analysis is originally expressed in terms of the absolute coordinates of samples. When the samples are distributed non-uniformly, or on a domain that is not a unit-hypercube, such analysis is difficult to carry. Differential analysis is instead expressed in terms of relative positions between two samples. By matching the distance function used to measure the distance between samples to the domain or the density of this distribution, it can produce a pertinent spectral analysis from non uniform sampling patterns or non Euclidean domains.

This differential formulation can be directly derived from the formulation of the Fourier transform for a set of samples. Formally, if we formulate the Discrete Fourier transform of a point set P n as follows

F ω = n-1 k=0 e -2πi(ω•x (k) ) (2.20)
The radial average is computed by averaging all the spectrum values between two s-D spheres, thus leading to a 1-D representation of the original s-D spectrum.

A pointset (right) with a finite number of frequency. Its power spectrum (middle) is clearly non isotropic (and was enhanced to make the peaks visible) but when you look at its radial average (left), the recurrent frquencies still appear clearly where x (k) is the k th sample from P n , ω is the frequency vector and • is the inner dot product. This leads to the following formulation for the Power Fourier Spectrum

P ω = |F ω | 2 = P r ω + P i ω (2.21)
where P r ω (resp. P i f ) is the real part (resp. the imaginary part) of P ω .

P r ω = 1 N n-1 k=0 cos(2πi(ω • x (k) )) 2
(2.22)

P i ω = 1 N n-1 k=0 sin(2πi(ω • x (k) )) 2 .
(2.23)

Using the sum rules of trigonometry, this leads to

P ω = 1 N n-1 j,k=0 cos(2πi(ω • (x (k) -x (j) ))) (2.24)
where P ω only depends on the differential between the pairs of samples x (k) and x (j) .

Therefore, using this type of analysis is just as discriminative as using traditional Fourier analysis. However, note that it applies efficiently only on domains where the distance function between pairs of samples can be properly defined for the entire domain (for example using geodesics for surfaces, of weighting the distance by the inverse density function for adaptive sampling).

This type of analysis prooves that the spatial organization of samples can be directly related to the spectral properties of this distribution. Furthermore, the results of spatial analysis can be strictly equivalent to the results of spectral analysis, when using proper tools, as presented in the next section.

Spatial Analysis

This spatial analysis of point sets is closely related to the spectral analysis. The simplest measure in spatial analysis is the minimal distance d (usually l 2 Euclidean distance) between two samples in a point set P n ,

d(x, x ) := min x,x ∈Pn ||x -x || 2 .
(2.25)

However, this measure, as simple as it is is not very useful. Instead, we will rely on the distribution of distances between pairs of points, called the pair correlation function.

Pair Correlation Function

In this whole section, we assume that our point set P n is stationary (see Section 1.4 on page 8). Before presenting the Pair Correlation Function (PCF), we will start by presenting two other second order estimators, known as the K-function and the L-function. For the notations, we will follow the choices made by [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF].

We will start by denoting N (i) (r), the function that estimates the number of samples in a radius r from a given sample x (i) (Fig. 2.10). This function in itself does not convey much information. However, it will be the basis for all second order spectral analysis functions, such as Ripley's K-function, the L-function or the Pair Correlation Function (PCF).

Ripley's K-function represents the average number of points that are at a distance d < r from a typical point. If we denote this value K(r), we have for a point set P n ,

K(r) := 1 λn n i=0 N (i) (r) (2.26)
where λ is the density of samples described in 1.4 on page 8. The normalization of K(r) by λ, allows K to be independent from the global density and thus to be more sensitive to local fluctuations.

This K-function relates the number of samples to the area of a ball of radius r. To alleviate this and relate instead the number of samples to r itself, one can use the L-function. This variant of the K-function is defined as

L(r) := s K(r) b s (2.27)
where s is the number of dimensions and b s is the volume on the unit sphere in R s . Both those functions can be seen as the cumulative distribution function (with some normalization) of the PCF. This PCF, noted (r), is considered as the best second order characteristic. Intuitively, it expresses the number of samples that are at a distance r for a typical sample (Fig. 2.10 on the preceding page). For r ≥ 0, it can be defined as (r) := K(r) db s r s-1 .

(2.28)

Oztireli [START_REF] Öztireli | Analysis and Synthesis of Point Distributions based on Pair Correlation[END_REF] devised a simplified estimator for this measure in the particular case of isotropic and stationary point processes. We therefore estimate the PCF of a pointset P n in the unit domain [0, 1) s using the following formula

(r) = 1 n 2 r s-1 i =j k σ (r -d(x (i) , x (j) )), (2.29)
where d(x (i) , x (j) ) is a distance measure between x (i) and x (j) . The factor k σ is used to smooth out the function. Oztireli relies on this smoothing to assume ergodicity for all sets. He uses the Gaussian function as a smoothing kernel, but one could use a box kernel or a triangle kernel instead. To compute a PCF, we use this estimator with 3 parameters, the minimal r, r min , the maximal r, r max and the smoothing value σ. Those values are usually chosen empirically. Note that as the number of samples increase, the distances between samples will be very different for similar distributions. To alleviate this, we normalize the distances in our estimations using the maximal possible radius for n samples ( [GM09], Eq (5) ). The results of all those estimators are illustrated Figure 2.11.

The issue when computing the values for , K and L, is that we have as input a point set defined in a bounded domain. However, those estimators implicitly assume an infinite set. This creates a bias since for points close to the border of the domain, their nearest neighbours outside the domain will not be taken into account. To avoid this, one can periodically tile or mirror a point set. One can also use more advanced estimators that take those effects into account. However, with a sufficiently big set and with sufficiently small radiuses, those edge effects can often be ignored.

The PCF is also a very interesting estimator as has been proven to be the spatial equivalent of the Fourier transform. You can find the mathematical relationship between both those measures in [START_REF] Kailkhura | Stair blue noise sampling[END_REF]. We recall that the differential analysis (Section 2.2.3 on page 23) was also proven to be equivalent to the Fourier Transform. This means that in most cases, the spatial and spectral analysis are interchangeable, and the use of one or the other is partly dictated by which will make the computations easier. However, the spectral analysis has been thoroughly studied and has been proven to be directly correlated to the variance in Monte Carlo estimators. This studies have been partly done for the spatial analysis but it is still a work in progress. The relationships between those characterizations (spatial analysis, spectral analysis but also uniformity) and the error in Monte Carlo integration will be detailed in the next section.

Relationship to error in integration

All the previously detailed characterisation was not randomly chosen. Those particular measures are of interest to us as they all relate to the error in the resulting MC/QMC estimator. In this section we see how we can use this relationship to devise sampling strategies that reduce this error as much as possible. Note that similarly to the discrepancy, for unbiased samplers, the integration error reduces as the number of samples increases. Each sampling pattern therefore has a convergence speed, that characterizes how the variance of the MC estimator using this sampling pattern changes as the number n of samples increases. Our reference worst case is the convergence speed of a pure random sampler, which is O(1/n) (Fig. 2.14 on page 32).

We remind here that we define the squared error in integration I e as Δ = |I n -I|, (2.30) where I is the real integration value and I n is the approximation obtained from estimating it with a MC/QMC estimator, using n samples obtained from a given point process.

Thanks to the Koksma-Hlawka theorem [START_REF] Hlawka | Funktionen von beschränkter variatiou in der theorie der gleichverteilung[END_REF], the uniformity of a point set happens to bound the integration error of QMC estimators. We detail this theorem first, and then, we see how in the case of stochastic patterns, the error can be reduced thanks to variance analysis. We remind that in this section, we focus on unbiased samplers, and therefore the error in integration is strictly equivalent to the variance (see Section 1.1.2 on page 4).

Koksma-Hlawka

The first theorem that interests us is the Koksma-Hlawka theorem [START_REF] Hlawka | Funktionen von beschränkter variatiou in der theorie der gleichverteilung[END_REF]. This theorem was developed by Hlawka and proven by Koksma (for the 1-D case) [START_REF] Jf Koksma | Een algemeene stelling uit de theorie der gelijkmatige verdeeling modulo 1[END_REF]. It bounds the integration error Δ using the discrepancy of a point set using the following inequality

Δ ≤ D * (P n )V (f ), (2.31)
where V (f ) is assumed to be finite, and represents the variation of the function f in the sense of Hardy and Krause (definition below) [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF].

This boundary implies that for point sets that are low discrepancy, the integration error is supposedly lower than for other point sets. An extension of this theorem was developed for the generalized l 2 discrepancy [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF] as

Δ ≤ D 2 (P n )V 2 (f ), (2.32)
where V 2 (f ) is finite and is the l 2 equivalent of V (f ).

As detailed in [START_REF] Keller | Myths of Computer Graphics[END_REF], this theorem is however a little unsuited for computer graphics. The issue is that it only applies on the class of functions of bounded variation in the sense of Hardy and Krause. A 1-D function is said to be of bounded variation if its total variation is finite. More formally, we denote V b a the variation of a 1-D function f in an interval [a, b] with

V b a (f ) := sup P ∈P n P -1 i=0 |f (x i+1 ) -f (x i )| (2.33)
which takes the supremum over the set

P = {P = {x 0 , x 1 , . . . x nP }| P is a partition of [a, b] with x i ≤ x i+1 for 0 ≤ i ≤ n P -1 } of all partitions of [a, b]. Then, f is a Bounded Variation (BV) function if V b a (f ) < +∞. ( 2

.34)

A BV function in the sense of Hardy and Krause is an extension of this definition to 2-D functions, which computes if the function is BV along the x and the y axis. In [START_REF] Keller | Myths of Computer Graphics[END_REF], Keller argues that this class of functions is very unhandy as, for discontinuities that are not axis aligned (which are very common in graphics), V b a becomes infinite. Following this observation, [Kel06] derived several new inequalities to relate the discrepancy to the integration error. However, as they imply diving into deeper details of BV functions and lead to results similar to the Koksma-Hlawka theorem, we will not detail this here. Interested readers can refer to [START_REF] Keller | Myths of Computer Graphics[END_REF].

We can already note that such theorem implies that a Low Discrepancy sampler, with a discrepancy of O(log(n) s-1 /n) when the number of samples n increases has a 2-D squared convergence speed of O(log(n) 2 /n 2 ) (Fig. 2.14 on page 32).

Even though those boundaries are very useful for deterministic samplers, they are still of limited use, since they are just an upper bound instead of a direct relation. This direct relation can be found when using stochastic samplers and variance analysis.

Integration error analysis

For stochastic samplers, we don't rely much on the integration error, but instead we rely on the MSE (as detailed Section 1.1.2):

MSE = Bias 2 + Variance.
(2.35)

Since we are only interested in unbiased stochastic samplers, this formulation states that the observed integration error is solely related to the variance of the approximation.

In the case where the error is solely due to the variance, we can rely on two very strong results. The first one [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF], directly relates the variance to the Fourier spectrum of both the sampling pattern and the function to sample. The second one [Özt16] relates the PCF of a sampling pattern to the integration error.

Spectral variance analysis

Many papers have conducted analysis of the variance in Monte Carlo integration, in both the Euclidean and spherical domains [Dur11; SK13; Hes+10; Ram+12; Pil+15]. We will detail here the results obtained by [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF]. Their result is an expression of the variance observed when using stationary point patterns, as a direct function of their Fourier spectrum. We denote I the function to integrate, and I n the MC estimation of the integral of I with a point set P n . The variance of I n , is given by

var(I n ) = μ(Ω) 2 n Ω P s (ω) P I (ω)dω (2.36)
where μ(Ω) is the Lesbegue measure for the unit toroidal domain [0, 1) s and P s , P I are the Power Spectra of respectively the sampling pattern and the integrand I. This formula can be further simplified for isotropic point patterns. Expressed in polar coordinates (n, ρ), we have

var(I n ) = μ(Ω) 2 μ(S s-1 ) n ∞ 0 ρ s-1 Ps (ρ) PI (ρ)dρ (2.37)
where S s-1 is the surface area of the s-D sphere, and Ps , PI are the radial average of the power spectra of respectively the point set and the integrand function. Note that as we are using the radial power spectrum, we assume an isotropic spectrum.

It has been shown that for natural images, PI will approximate a 1 ω curve (see Chapter 7). This means that as ρ increases in Equation 2.37, the value of PI will tend to 0. This means that to reduce the variance, the value of Ps must be as small as possible when ρ is small (Fig. 2.12 on the facing page).

Our goal is therefore, for integration purposes, to use preferably sampling patterns that present as little energy in the low frequencies as possible. In 1983, Yellot [START_REF] John I Yellott | Spectral consequences of photoreceptor sampling in the rhesus retina[END_REF] observed that the photoreceptors in the retinas of monkey where distributed following a very uniform but yet irregular distribution. It has been theorized by Ulichney [Uli88] that this distribution, called Blue Noise, leads to the the best spectrum for integration. This ideal spectrum presents no energy in the lowest frequencies, followed by a peak and a flat region (Fig. 2.13). Currently, samplers that generate samples with such a spectrum are referred to as Blue Noise samplers.

Following this analysis, the theoretical MSE of Monte Carlo estimators using a stratified sampler decreases in O(1/(n √ n)) as n increases. A similar convergence speed was also empirically devised for the Blue Noise samplers. This is better than the speed of the pure random sampler but asymptotically not as good as the convergence speed of Low Discrepancy Samplers (Fig. 2.14 on the following page).

Spatial variance analysis

In [Özt16], the authors devised a formula relating the PCF of a stationary sampling pattern to its resulting variance in MC integration. We will denote the PCF of a sampling pattern. From Campbell's theorem, they expressed the expected value of the MC integration of a function f defined over a domain Ω using the statistics of the sampler used. If we denote I n the MC estimation of the integration of f over Ω, this leads them, for stationary processes, to The convergence speeds of the main categories of samplers, when the number n of samples increases. In 2-D, the pure random sampler, which is our worst case reference, converges in O(1/n). The samplers with as little energy as possible in their low frequencies converge in O(1/(n √ n)), when the low discrepancy samplers converge in O(log(n) 2 /n 2 ). where λ is the intensity function (constant for stationary processes). They also derive E(I n ) 2 as

E(I n ) = Ω f (x)λdx (2.
E(I n ) 2 = Ω f 2 (x)λdx + E( i =j f i f j ) (2.39) = Ω f 2 (x)λdx + Ω×Ω I(x)I(y)λ 2 (x -y)dxdy. (2.40) Finally, from var(I n ) = E(I n 2 ) -(E(I n )) 2 (2.41)
they derive the formula relating the variance to the PCF,

var( În ) = 1 λ f 2 (x)dx + a f (h) (h)dh - f (x)dx 2 (2.42)
where h = xy, and a f represents the auto-correlation of f , which tends to 0 as h increases due to the finite support of f .

If we look at each term of this function, we see that the first term depends on 1 λ , meaning that the more samples the higher the intensity and the lower the error. The third term does not involve the point pattern, so we only care for the second term. What can be seen here is that when h is small (when a f is not 0), the value (h) should be as small as possible. This underlines the spectral analysis from [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF] that stated that as the Power spectrum of the function tends to 0 as the frequency increases, the Power spectrum of the point set should be as small as possible in the low frequencies. It also underlines the correlation between the PCF of a set and its radial spectrum.

Conclusion

In this section, we described several measures used to characterize a point set. Among them are the measures describing the uniformity of a set. The main one is the l ∞ star discrepancy. Thanks to the Koksma-Hlawka theorem, we know that a point set with a lower discrepancy limits the maximal possible error. As this value can be difficult to compute in s-D, one can rely on the generalized l 2 star discrepancy to ease the computations while preserving the Koksma-Hlawka inequality.

Then, we have all the measures that rely on computing statistics over the spatial distribution of samples. Those measures have been shown to be roughly equivalent to the spectral properties of the point set, while being in any cases much easier to compute, even in s-D. Furthermore, they are also directly related to the integration error, which makes them our main criteria when evaluating sampling patterns.

In the following chapter, we use all those characteristics to classify existing sampling techniques and sampling patterns, and to derive their convergence speeds.

http://cpc.cx/jDN http://cpc.cx/jEW 3

State of the art

In this chapter, we present existing samplers, evaluating them based on the characteristics detailed in the previous section. We classify those samplers into 3 categories. First, in Section 3.1 on the next page, we present some historical samplers. Such samplers usually aim at improving over the pure random sampler w.r.t. its uniformity. They are usually very simple and do not provide strong guarantees over the spatial organization of samples. However, such samplers are very efficient, and extend easily to higher dimensions. This means that even though, they are not optimal for Monte Carlo Rendering (where a higher correlation between samples is required to ensure a uniform and unstructured coverage of the integration domain), they are still implemented in most current rendering systems [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF]. Section 3.2 on page 40 presents the samplers that generate Poisson Disk or Blue Noise spectrum. Both those types of samplers increase the amount of correlation between samples to try generating stochastic sets of samples with roughly the same distance between each sample and all of their neighbours. Note that the Poisson disk Spectrum differs slightly from the Blue Noise spectrum, as it presents a non-zero energy in the lowest frequencies, which is harmful to the variance of MC estimators (Section 2.4 on page 28). Still, in a vast part of the literature both those spectra are considered the same which is why we regrouped them. Contrary to historical samplers, Blue Noise and Poisson Disk samplers are limited to low dimensions (usually, they only generate 2-D samples).

Finally, Section 3.3 on page 57 details the samplers that are low discrepancy. Such samplers are usually defined from implicit analytical functions, based on bit reversal. This makes them very efficient, even in high dimensions, and allows to add consistent new samples iteratively within an existing point set. However, they are rarely stochastic and their analytical definition generates highly structured sampling patterns. Note that in most of this section, the samplers will be defined over the unit domain [0, 1) s (unless specified otherwise). Figure 3.37 on page 71 presents a summary of the main samplers from each category, evaluated on the closeness of its spectrum to Blue Noise and on the quality of its discrepancy. Figure 3.38 on page 72 presents the same samplers, evaluated on their computing efficiency (speed) and on the number of dimensions that can be sampled.

Historical samplers

In this section, we present some historical samplers. Most of those samplers aim at improving over the pure random sampler. Thus, they usually present a higher error in integration than optimal Blue Noise samplers, but they do not generate aliasing during rendering since they do not present any structures in the sampling patterns. They are also very easy to implement and efficient even in higher dimensions. This is why they keep being so popular in rendering engines [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF].

Whitenoise sampler We start by presenting the most naive sampler, the Whitenoise sampler. This sampler generates purely random samples. Consequently, in the resulting point set, all frequencies are present with the same energy which leads to a flat Fourier Spectrum (Fig. 3.1). Note that since the Fourier spectrum is flat, there is a higher energy in low frequencies, which leads to a high variance in integration. Due to the simple expression of this sampler, the MSE of the pure random sampler is known to be in O 1 n , with a discrepancy in O 1 √ n .

Stratified sampler One of the drawbacks of using a Whitenoise sampler is that, as it is fully random, nothing prevents having subsets of samples that are arbitrarily close to each other. This translates as portions of space that contains a lot of samples, and others that are almost empty (see Figure 3.2 on the facing page, left). As the energy in the low frequencies of the spectrum is due to those empty spaces, we want to avoid them as much as possible. One way to do that, is to use Stratified sampling. This sampler, in its simplest form, partitions the space using a s-D grid, and generates a random sample inside each cell of this grid (see Figure 3.2 on the next page, right). This limits the number of arbitrarily close samples to subsets of size 2 s and thus increases the uniformity. As can be seen in the radial spectrum (Fig. 3.3 on the facing page), there are much less energy in the low frequencies which leads to a lower variance in MC integration, as detailed in [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF] (Section 1.1.2 on page 4).

This sampler is one of the most used in rendering engines as its variance is reasonably low and its implementation is very easy. However, this sampler cannot generate any number of samples. In dimension s, this sampler can only generate a point set P n when n = k s . Therefore, as we want to integrate over more and more dimensions, we will have to generate more and more samples just to be able to properly sample each dimension. This phenomenon is called the curse of dimensionality, and is the main drawback of the Stratified sampler relatively to the Whitenoise sampler. This is usually alleviated in rendering by generating 2-D stratified sampling patterns that are then shuffled between dimensions [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF].

Whitenoise Stratified

This sampler has an improved uniformity but, several samples can still present the same x or y coordinate, leading to overlapping samples in the 1-D projections. As detailed in Section 2.1 on page 14, this overlap can increase drastically the discrepancy of a point set and thus should be avoided. To alleviate this, sampling strategies based on shuffling a canonical point set have been developed.

Shuffling canonical projections Most those methods aim at generating stochastic point sets while preserving the (s -1)-D projections [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF][START_REF] Chiu | Multi-jittered sampling[END_REF][START_REF] Kensler | Correlated multi-jittered sampling[END_REF]. As explained in Section 2.1 on page 14, having two samples on the same line or on the same column is very harmful to the uniformity of this pattern. In the context of the samplers described here, preserving the (s -1)-D projections of a point set is used to ensure that those projections will be stratified. This supposedly decreases the error in the Monte Carlo estimators and therefore improves the quality of the rendered image. This property is illustrated Figure 3.4.

Those samplers usually start with a canonical set, either stratified or diagonalized (Fig. 3.5 on the next page). Then, they randomly swap the lines and columns in this arrangement. This preserves the 1-D projections, as no samples were on the same line or column originally and the coordinates are not changed, merely swapped.

The N-Rooks sampler is one of the main examples of such samplers [START_REF] Shirley | Discrepancy as a quality measure for sample distributions[END_REF]. But even though it preserves the stratification of the 1-D projections, it presents spectral characteristics that are almost identical to Whitenoise (Fig. 3.6). The only notable change is the presence of a "black cross" in the middle of the Fourier Spectrum (since there are no samples aligned on either dimension, there are no frequencies in the strictly vertical or strictly horizontal directions). The original goal of this sampling was to allow a stratified sampling for any number of samples, instead of simply for point set whose size is a power of s, but due to the poorer quality of its spectrum, this sampler is almost unused in rendering nowadays.

One can note that some authors [START_REF] Kensler | Correlated multi-jittered sampling[END_REF] alleviate this "Whitenoise effect" by correlating the swapping of coordinates. Their goal was to combine the very good properties of the stratified sampling with a higher uniformity, ensured by their original canonical arrangement. However, even if this gives good results for low number of samples, the correlation generates spectral peaks for a higher number of samples (Fig. 3.7 on the next page).

If we look at all the samplers presented in this section, we can note that are all very efficient and easy to implement, even in higher dimensions. They are therefore quite present in current rendering systems. However, as they rely or pure randomness (either in the sample generation or in the shuffling), they do not present Fourier spectra that are as good as the spectra from Blue Noise sampling patterns (Section 3.2 on the following page), or the excellent uniformity of Low discrepancy samplers (Section 3.3 on page 57).

Poisson Disk and Blue Noise samplers

Using the analogy with the visible light spectrum, several colors of noise were devised (Fig. 3.8 on the facing page). The most famous one is the white noise, generated from pure randomness with a flat spectrum. One can also note the existence of red or pink noise, that present more energy in the lowest frequencies. Such noise are formed by clustered points and are unsuited for Monte Carlo integration, but are very useful in biology or physics to represent the distributions of for example plants or galaxies. Green noise presents higher mid frequencies, and can be used to modelize fallen leaves or the distributions of seeds around trees for example. Finally, Blue Noise puts the emphasis onto high frequencies. This type of noise is of major interest for Monte Carlo integration and is the kind of distribution that interest us here. A good survey of Blue Noise sampling techniques can be found in [START_REF] Yan | A survey of blue-noise sampling and its applications[END_REF]. Many notations in this section will come from this survey.

Poisson Disk samplers

The most naive sampler generating Poisson Disk distribution is the Dart Throwing sampler [START_REF] Robert | Stochastic sampling in computer graphics[END_REF] (Fig. 3.9 on the facing page). The idea is to iteratively insert samples into an already existing set. A new sample is generated purely at random, and is kept only if it is at a sufficient

White Noise

Pink Noise

Green Noise Blue Noise distance from every other samples. This distance is called the disk radius, and is computed based on the maximal packing density γ s of n s-D spheres with the following formula

r max := s γ s n Γ( n 2 + 1) π n/2 , (3.1)
where the Γ function is the extension of the factorial function to positive complex and real-valued arguments, defined as

Γ(x) := ∞ 0 x z-1 e -x dx (3.2) with x ∈ C.
However, this algorithm has a major flaw: it is not guaranteed to end. If one wants to generate n samples, it is possible to reach a situation where m < n samples where placed, but where there is no longer any position in space that is at a sufficient distance from all other samples. To alleviate that, the relaxation dart throwing method [START_REF] Mccool | Hierarchical Poisson disk sampling distributions[END_REF] was developed (Fig. 3.10). This algorithm initiates the radius r max with a large value, and if there is too many unsuccessful trials, this radius is decreased. This process goes on until all n samples have been placed. Note that this relaxed dart throwing can be considered as a sequential sampler, as it iteratively adds samples within a pre-existing set, while being consistent regarding the distance between pairs of samples.

Note that those two methods, generate different types of Poisson disk sampling. The dart throwing sampler belongs to the family of accurate samplers, which generates point sets such that all the distances d between sample pairs are d > r max . In opposition, the relaxed dart throwing algorithm is an approximate sampler. It generates point sets where the distance between the majority of sample pairs is d > r max but where there might be pairs of samples with d < r max . Another side effect of the relaxed Dart Throwing is that the resulting sampling pattern is not necessarily maximal. This means that samples could still be added to the set without changing the value of r max .

For both those algorithms, one can note that they are highly inefficient. You may generate lots of samples before finding one that fits. Therefore, another simple approximate algorithm for Poisson disk sampling was devised, the Best candidate algorithm [START_REF] Mitchell | Spectrally optimal sampling for distribution ray tracing[END_REF]. This algorithm generates a new random samples by generating k new candidates, and then by selecting the sample with the largest minimum distance to all others. Note that this algorithm is guaranteed to finish and is much faster than dart throwing. However, it is sill has an O(n 2 ) complexity.

Another noticeable approach, although quite exotic, is to generate Poisson Disk through molecular dynamics [START_REF] Boris | Geometric properties of random disk packings[END_REF]. Such methods generates n random particles and moves them around with a collision check. The radius of the particles gradually increases until no particle can move anymore. The main asset of such methods is that particular dynamics offers well known and efficient algorithms, that can be parallelized on GPU. This makes such approach faster than the previous ones, but it may generate very regular structures in 2-D. Other similar methods will be devised for Blue Noise at the end of Section 3.2.2 on page 44. We also point that [START_REF] Wei | Parallel Poisson disk sampling[END_REF] can generate higher dimension poisson disk sampling, although its practicability decreases as the dimension increases.

In his paper, Bridson [START_REF] Bridson | Fast Poisson disk sampling in arbitrary dimensions[END_REF] accelerates the generation of samples in high dimension by partitioning the domain using a regular grid. This grid resolution is carefully chosen so that each cell contains at most one sample. Then, when generating a new sample, we use the best candidate algorithm but only test the distance against the samples in the neighbouring cells instead of testing the distances against all samples. This algorithm runs in O(n) but generates approximate distributions. There exists similar algorithms that runs in O(n log(n)) for accurate distributions [START_REF] Manuel | Accurate multidimensional Poisson-disk sampling[END_REF]. However, since all those methods rely on a grid subdivision of the domain, they all face the curse of dimensionality. This means that as the number of dimension increases, the number of grid cells increases exponentially. To alleviate that, [START_REF] Ebeida | A Simple Algorithm for Maximal Poisson-Disk Sampling in High Dimensions[END_REF] uses a flat implicit representation of a quadtree that is also ported on the GPU.

Note that until now, most of the aforementioned Blue Noise methods aimed at generating 2-D (or low-D) Poisson Disk sampling patterns. This limitation comes from the sampling structures used that were defined in 2-D, or simply to the computational cost of the algorithm that prevents its use in higher dimensions. For all those methods, we still have a non-zero energy in the low frequencies (Figures 3.10 on page 42, 3.11 on the previous page), in [START_REF] Torquato | Random sequential addition of hard spheres in high Euclidean dimensions[END_REF], authors found that, for a point set P n with sufficiently large n, the minimal value of the Power Spectrum near the DC was 1/20. It will therefore have a very bad behaviour regarding variance reduction, due to the relationship between the spectrum and the variance [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF] (see Section 2.4.2 on page 30). Figure 3.12 illustrates the product of the radial spectrum of a Poisson Disk distribution [START_REF] Dunbar | A spatial data structure for fast Poisson-disk sample generation[END_REF] with the radial spectrum of a typical scene. One can note that since the radial spectrum of the point set is not null in the lower frequencies, the product has a higher value near 0. Since the variance is directly linked to the integrand of this product [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF], this small peak at 0 automatically leads to a higher variance. Poisson Disk are therefore mainly used for stippling or object placement but not for Monte Carlo integration, which is our topic of interest. Instead, we prefer to use Blue Noise samplers.

There are several families of methods to generate Blue Noise samplers. The main ones are the algorithms based on energy minimization (Section 3.2.2), and the tiled based systems (Section 3.2.3 on page 49). But Blue Noise can also be generated from dithering matrices (Section 3.2.4 on page 52).

Optimization based Blue Noise

The samplers discussed in this section present the ideal spectrum for variance reduction as they present no energy in the lowest frequencies (Fig. 3.8 on page 41). We refer to this energy-less interval as the zero region of a spectrum. In this section, we present methods that generate Blue Noise from energy minimization. Note that most Blue Noise samplers were initially devised to account for adaptive sampling. We remind the reader that in this document, we only focus on uniform sampling, but we will still compare the performances of Blue Noise samplers regarding their adaptivity. As it was one of their original purpose, this whole section would make little sense if we did not. Nonetheless, our main focus will stay on their variance reduction property.

The simplest algorithm to generate Blue Noise from iterative optimization is known as Lloyd's relaxation [START_REF] Lloyd | Least squares quantization in PCM[END_REF][START_REF] Du | Centroidal Voronoi tessellations: Applications and algorithms[END_REF], also referred to as Centroid Voronoi Tessellation (CVT) in Computer Graphics (Fig. 3.13). This algorithm generates the Voronoi diagram of a set of samples and displace each sample to the center of its Voronoi cell. The process is then repeated, until the system reaches stability (Fig. 3.14). Using such methods, a point set P n is optimized by minimizing the following energy

E CV T := n i=1 Vi ||x -x i || 2 dx (3.3)
where {V i } n i=1 is the set of cells of the Voronoi diagram with

V i := {x ∈ Ω, ∀j||x -x i || ≤ ||x -x j ||}. ( 3.4)
A key to the efficiency of such algorithms is the computation and update of the Voronoi diagram. In this document, we will assume that this computation is done and we will not go into further details regarding how it is done. An interested reader can refer to [All+05; Val+08; Yan+09; Ron+11; Yan+14b] for more information regarding Voronoi diagrams. In 2-D, point sets computed using Lloyd's algorithm converge to an hexagonal grid. To alleviate this, the optimization must be stopped after a certain number of iteration, and before reaching convergence. Unfortunately, there is no objective rule as to how to choose this number of iterations [START_REF] Balzer | Capacity-constrained point distributions: a variant of Lloyd's method[END_REF]. Furthermore, the computation of the Voronoi diagram limits this algorithm to low dimensions.

Early algorithms aimed at generating Blue Noise samples for stippling applications. Therefore, their goal was not the Fourier spectrum of the resulting set but its visual quality and adaptivity. Their method was then improved upon by [START_REF] Li | Fast capacity constrained Voronoi tessellation[END_REF] that proposed an accelerated algorithm to construct capacity constrained Voronoi tessellations, allowing this computation to be parallelized on GPU. Another paper by Xu et al. [START_REF] Xu | Capacity-constrained Delaunay triangulation for point distributions[END_REF] aimed at speeding up this Blue Noise generation. In this paper, the authors no longer use a CCVT but instead the Capacity Constraint Delaunay Triangulation (CCDT), as the Delaunay triangulation is the dual structure from the Voronoi diagram. This method is faster than [START_REF] Balzer | Capacity-constrained point distributions: a variant of Lloyd's method[END_REF] and extends the definitions to mesh surfaces. Finally, [START_REF] Chen | Variational blue noise sampling[END_REF] combined CCVT and CVT for sample surfaces. Their energy function is

E CapCV T := E CV T + λE CapV T (3.7)
where E CapV T is defined as

E CapV T := i=1 n Vi ρ(x)dx . (3.8)
They also achieve a significant performance improvement from an efficient optimization framework.

Other methods were devised based on [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF]. This method places samples one by one, placing each new sample in the region with the less samples already present. They define the position of this sample as the point of space that is further from any other samples. Proof exists that this sample necessarily lies on the Voronoi diagram of the current set. Therefore, at each iteration, the algorithm tests all points on the edges of the Voronoi diagram of the set, adds a point on the one that is farthest from all samples, and updates the Voronoi with this new sample. This idea was extended in [START_REF] Schlömer | Farthest-point optimized point sets with maximized minimum distance[END_REF], with their Farthest Point Optimization (FPO) method. In this paper, the whole set of samples is generated at random and iteratively optimized. Each iteration selects a sample, and pushes it as far away as possible from all other samples. However, the resulting pattern still presents some regularity that we would like to avoid (Fig. 3.16). In [START_REF] Kanamori | Deterministic blue noise sampling by solving largest empty circle problems[END_REF] the authors presented a similar method using a Delaunay triangulation updated iteratively. Using this Delaunay triangulation, the FPO algorithm was parallelized on the GPU by [START_REF] Chen | Parallel Blue-noise Sampling by Constrained Farthest Point Optimization[END_REF]. Finally, it was extended to surfaces and adaptive sampling by [START_REF] Yan | Blue-Noise Remeshing with Farthest Point Optimization[END_REF].

A very notable method for Blue Noise generation was developed using Optimal Transport (BNOT) [START_REF] De Goes | Blue noise through optimal transport[END_REF] (Fig. 3.17 on the following page). They use a Power diagram instead of a Voronoi diagram. A Power diagram is an extension of the Voronoi diagram where each cell is defined as

V w i := {x ∈ Ω, ||x -x i || 2 -w i ≤ ||x -x j || 2 -w j ∀j} (3.9)
with W = {w i } n i=1 the set of weights defined for each sample. They then minimize the following function

E BN OT = (W ) + i=1 nλ i (c i -c) (3.10)
where

c i = V w i ρ(x)dx) and (W ) is defined as (W ) = i=1 n V w i ρ(x)||x -x i || 2 . (3.11)
There also exist methods that generate Blue Noise with particle based approaches [START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF][START_REF] Jiang | Blue noise sampling using an SPH-based method[END_REF]. In [START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF] (Fig. 3.19 on page 49), the authors motivation is to have a better computational efficiency than optimization based methods. They attach a density function to each sample and measure the difference between the sum of densities from the samples and the ex- pected density. This gives an energy function, but instead of minimizing this energy, they use it to define a statistical model, from which they draw points. This kind of algorithm is slow, but by generating samples one by one, the authors manage to generate their samples in linear time.

It is worth noting that this algorithm was initially thought for stippling, and cannot be used for Monte Carlo integration. The issue is that, as the authors draw samples from statistical models, there is a minor shifting of the samples, imperceivable for stippling, but that generates a lot a bias when used with Monte Carlo estimators.

Other particle based approach are the ones of [START_REF] Jiang | Blue noise sampling using an SPH-based method[END_REF]. Their general algorithm is illustrated Figure 3.18 on the facing page. The authors rely on Smoothed Particle Hydrodynamics (SPH), a method from fluid simulation that minimizes the internal pressure variance. This method uses a spring mass system if the difference between the density of particle and the rest density is positive, the particles are pushed away, and if the difference is negative, the point are brought together. The adaptation of [START_REF] Jiang | Blue noise sampling using an SPH-based method[END_REF] is to have clamped this difference to 0, so that the points are never brought together. Each sample is associated to a density kernel. At each iteration a sample is added and the spring mass system is updated. Their method allows to generate higher dimensions Blue Noise at faster rates than other optimization based methods. More recently, [START_REF] Abdalla | A simple push-pull algorithm for blue-noise sampling[END_REF] proposed a new particle based Blue Noise method, aiming at optimizing the coverage radius while preserving computational efficiency. However, to the best of our knowledge, such methods usually focus on surface sampling for remeshing and the resulting sets where never thoroughly tested for Monte Carlo integration.

All aforementioned methods tend to generate Blue Noise in "low-D". However, in a rendering context, we need to solve several recursive 2-D integrands (see Section 1.1.1 on page 2). In this context, the method in [START_REF] Reinert | Projective Blue-Noise Sampling[END_REF] is notable as its aim is to generate s-D sampling patterns, with the added constraint that they present the same spectrum in every projections. Authors suggest two methods. The first one, is derived from the dart throwing algorithm, which generates Poisson Disk sampling patterns in every projections. To do so, they add new constraints on the Dart Throwing algorithm, one for each projection. As in the original Dart Throwing algorithm, a new sample is added if it is a proper distance from all other samples in the domain, but it is also tested over all samples within each projections of the s-D set. Only a sample passing all those tests is added to the final set. The second method relies roughly on the same idea of adding projection-wise constraints, but is derived from the Lloyd's relaxation algorithm, and generates Blue Noise sampling patterns in the projections. Unfortunately, despite the s-D claim of the paper, both those algorithms are so costly that they hardly generate samples in more than 4 or 5 dimensions. They are also almost unable to generate a high number of samples.

Most their examples are limited to a few thousands samples.

All the methods in this section rely on iterative optimisations, usually through the computation and update of Voronoi Diagrams or through particle simulation. Those optimizations generate high quality Blue Noise but it has a computational cost. To alleviate this cost, Tiled Based methods were devised, based on precomputed tiling patterns, that simply needs to be combined together at run-time.

Tiled Based algorithms

In this section, we focus on methods generating Blue Noise from precomputed tiles filled with samples. Such methods are very efficient as they are based on simple recursive paving rules and tiles combinations. The main difficulty is avoiding artefacts due to the repetitive/periodic

Set of 11 wang Tiles

Paving from those tiles Wang tiles can be used for an aperiodic paving of space (meaning without periodic structure) with square tiles (note that those tilings are not mathematically aperiodic, but as the tiling is randomly constructed, there is a probability 1 that it will not be periodic). To handle the boundary issues, each edge of each square tile is associated with a color. The final tiling can then only be composed of tiles with similar facing edge color, under the constraint that no tile can be rotated.

To use this tiling for Blue Noise, one fills the tiles with precomputed sets of Blue Noise samples. Those sets can be generated by a constrained Lloyd's relaxation, where this relaxation is additionally constrained on the edges to maintain a coherency between edges of the same color. Amongst such methods, the most notable is [START_REF] Kopf | Recursive Wang tiles for real-time blue noise[END_REF]. This algorithm, contrary to [Hil+01; Coh+03; LD05], supports adaptive sampling and random direct access of samples. Furthermore, the precomputation of their tiles takes only a few minutes and weight only a few mega bytes of storage. In their experiments, authors claim to generate a set of 8 tiles with 2048 samples each in around 20 minutes, stored in a 6 megabyte file.

The drawback of this method is that it uses several samples per tiles. This implies that to generate an arbitrary number of samples, one needs to rely on a ranking algorithm. This ranking [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF] is not an off-line process and does not allow a fine control for adaptivity. Furthermore, storing the tiles or the recursive paving rules can be memory costly or they may be too little available tiles to avoid repetition in Monte Carlo rendering context.

Other methods use more complex tilings with a single sample per tile (Fig. 3.21 on the next page). This allows for a finer control for adaptivity and avoids the memory cost of tile storage. However, such methods are even more prone to aliasing than the previous ones. This is the case of Penrose tiling [START_REF] Ostromoukhov | Fast hierarchical importance sampling with blue noise properties[END_REF]. This method presents a recursive definition of a Penrose tiling where each tiles contains a sample. Furthermore, a set of precomputed correctional offsets (from a Lloyd's relaxation) is used to ensure Blue Noise properties for the final set. However, the obtained set presents several peaks in its Fourier spectrum (Fig. 3.22 on the facing page).

To alleviate that, one can use polyominoes to tile the domain [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]. The process of sampling with polyominoes is almost identical as sampling with penrose tiles, with a different correctional

Penrose Tiling [Ost+04]

Polyominoes Tiling [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF] Irregular Trihexes tiling [Wac+14] offsets table. Even though sampling with polyominoes effectively reduces the amount of aliasing, it does not remove it entirely.

Through the usage of irregular trihexes, [START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF] achieved a tiling that is both recursively subdividable and without any obvious regularity. The resulting sampling set thus have a Fourier spectrum that is without peaks. Furthermore, their method also allows for user defined Fourier spectrum. The Fourier spectrum is controlled by the set of correctional offsets for the samples positions in the tiles. By changing such offsets, the final Fourier spectrum can be changed. However, due to its irregularity, their tiling patterns requires a table storing the subdivision rules. Combined with the table containing the offsets, this algorithm becomes very memory costly. All together, the tables given by the authors for the subdivision rules and the Blue Noise offsets weight about 2 gigabytes.

Finally, a notable method for Tiled based Blue Noise is the use of AA Patterns [START_REF] Abdalla | AA patterns for point sets with controlled spectral properties[END_REF] (Fig. 3.23 on the next page). As noted by [START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF], the tile storage can grow up to gigabytes of memory. Here, the authors aim at using tiles with an implicit definition, AA patterns. Those patterns come from error quantization and are used for ornamental purposes. By using those patterns as tiles and storing only the displacement vectors within each tile, AA patterns allow for tiled based Blue Noise generation without any memory costs. Their spectrum is also peak free due to the aperiodic tiling resulting from AA patterns. However, they no longer allow for adaptive sampling. Very recently, [START_REF] Ahmed | An Adaptive Point Sampler on a Regular Lattice[END_REF] presented a method that generates Blue Noise from a regular grid. This scrambling allows for a very efficient handling of the recursivity and of neighbouring tiles.

Thanks to their recursive paving, all tiled-based methods are naturally suited for adaptive sampling. An important number of those methods where originally calibrated for stippling purposes. This means they aim at representing an image a set of dots, which can have applications, for example, for printing purposes. Another class of algorithms that can be used for the same purpose, the dithering methods. Such methods will be the topic of the following section.

Dithering based methods

Originally, dithering was used when printing images. Printers only place ink dots on paper, therefore, a black and white image must be stippled (halftoned) to be printed. Dithering is used to prevent regularities in the halftoned image. There are several dithering algorithms [START_REF] Ulichney | Digital halftoning[END_REF]. The easiest one is to simply threshold the original image, meaning that ink is placed for every point that is above a given threshold (Figure 3.24 on the facing page (a)). However, with such method, there is a lot of detail loss. To avoid that, this threshold can be made random, but this results in a lot of noise (Figure 3.24 on the next page (b)). Instead, one can use ordered dithering. A dithering matrix is precomputed and used as threshold (Figure 3.24 on the facing page (c)). Among such matrices, some are tuned for Blue Noise sampling patterns. This is such matrices that will interest us here (Figure 3.24 on the next page (d)).

Among the first methods to do Blue Noise dithering is [START_REF] Mitsa | Digital halftoning technique using a blue-noise mask[END_REF]. They iteratively build dithering matrices from the difference between the Fourier transforms of the current matrix and the new matrix. Following this work, [START_REF] Purgathofer | Forced random dithering: improved threshold matrices for ordered dithering[END_REF] tried to generate such matrices with an algorithm very similar to dart throwing. Still another method was proposed by [START_REF] Jeffrey | Generation of blue noise arrays by genetic algorithm[END_REF] to generate such matrices from genetic algorithms. Here, the method that interests us more is the very recent method from [START_REF] Cornel | Forced Random Sampling: fast generation of importance-guided blue-noise samples[END_REF] (Fig. 3.25 on the facing page). They pre-compute matrices using the method from [START_REF] Purgathofer | Forced random dithering: improved threshold matrices for ordered dithering[END_REF] and then tile those matrices to generate a point set of arbitrary size. .

With such dithering algorithms, we completed our tour of the main families of methods generating Blue Noise sampling patterns. However, this list is in no way exhaustive (far from it).

Other notable algorithms exists that did not directly aimed at generating Blue Noise sampling distributions but that still may be considered as Blue Noise samplers. In the next section, we will detail some of those algorithms.

Other methods

In this section, we present methods that generate variants of the classical Blue Noise sampling patterns. Some try to generate anisotropic or Multi-class Blue Noise, when other aim at targeting a user defined Fourier spectrum, or improving over the original Blue Noise spctrum.

Anisotropic and Multi Class Blue Noise

Those methods are of little use for Monte Carlo integration. However, they are notable variants of the Blue Noise sampling patterns and are presented here for the sake of exhaustivity. The anisotropic Blue Noise presented in [START_REF] Li | Anisotropic blue noise sampling[END_REF], is generated from an adaptation of the Dart throwing and the Lloyd's relaxation algorithm. The idea is to generate samples along a given vector field, mostly for stippling and importance sampling of surfaces (Fig. 3.26). For object placement, in many situations you may have to organize several families of objects. In [START_REF] Wei | Multi-class blue noise sampling[END_REF], the authors present a method for Multi-class Blue Noise. They also adapt the Dart throwing and Lloyd's relaxation to handle the attachment of an identifier to each sample, identifying the class of this sample.

Controlling the Fourier spectrum

A few methods were devised allowing to generate samples matching a user-given Fourier Spectrum, usually based on methods like Dart throwing or Gradient descent optimization. The versatility of such algorithms is very important as it could help creating tailored sampling patterns when rendering a given known scene.

The first paper to note is the one from [START_REF] Öztireli | Analysis and Synthesis of Point Distributions based on Pair Correlation[END_REF] This idea of controlled Fourier spectrum was also developed in [START_REF] Zhou | Point sampling with general noise spectrum[END_REF], where instead of using PCFs, they use differential analysis from [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF] and a gradient descent analysis. Similarly to [START_REF] Öztireli | Analysis and Synthesis of Point Distributions based on Pair Correlation[END_REF], as their fitting is based on a isotropic measure, they find it difficult to generate structured sampling patterns.

Finally, we can note the recent paper of [START_REF] Roveri | General Point Sampling with Adaptive Density and Correlations[END_REF], that uses the PCFs to generate sampling pattern with adaptive density and correlations. Authors provide a framework allowing to perform consistent analysis of such sets. However, even though such adaptive correlation can be very useful for stippling, object placement or reconstruction, it is still unknown how they behave on a context of Monte Carlo rendering.

We point out that all tiled based methods allow for a certain amount of spectrum control by changing the offset table.

Variants on the Blue Noise spectrum

Other sampling patterns exist that aim at a particular Fourier spectrum, improving over the original Blue Noise spectrum. However, for the following methods, the aim is no longer to allow user control of the Fourier spectrum but to generate a point set with a Fourier spectrum that was chosen to have the largest possible zero region or to avoid the oscillation in the higher frequencies (or both). The first of such methods is the Step blue noise [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF]. The authors aimed at obtaining a radial spectrum in the shape of a step function (Fig. 3.27). Here, the authors mostly aim at having a flat final spectrum. This is motivated by the idea that the high amount of energy in some frequencies and the accompanying oscillations are the part of the spectrum that generate most of the noise in the resulting image. In the same paper, they also propose the Single-Peak blue noise function, that instead of presenting a range of frequencies with higher energy, presents a single higher peak (Fig. 3.28). In both cases, the spectra are obtained from a gradient optimization.

A very recent paper [START_REF] Kailkhura | Stair blue noise sampling[END_REF] tried to increase the size of the zero region to its theoretical maximal size. To do so, they developed the Stair Blue Noise sampling pattern. They generate such a set from a target PCF, and use a least square optimization to optimize the PCF of a given set to match the Step Blue Noise target. The asset of this method is that it allows for a larger zero region than other methods, which reduces the variance in Monte Carlo integration [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF].

We recall that in this document, we are focusing on unbiased samplers, therefore, in the context of this PhD, reducing this variance is equivalent to reducing the integration error. It is this property that motivates research for more efficient and higher dimensional Blue Noise sampling patterns. However, there are other ways to reduce this integration error. Another possibility is to rely on sampling patterns with a mathematically optimal uniformity. Such samplers are referred to as Low discrepancy (LD) samplers, and will be the topic of the next section.

Low discrepancy samplers

A sampler is usually said to be Low discrepancy if its l ∞ star discrepancy in s dimensions decreases in O(log(n) s-1 /n) when the number of samples n increases. Most of those sets are deterministic so it makes no sense to evaluate their variance properties. Instead, we simply evaluate their integration error, using the theorem of Koksma-Hlawka [START_REF] Hlawka | Funktionen von beschränkter variatiou in der theorie der gleichverteilung[END_REF]. This theorem states that the discrepancy of a sampler is an upper bound over its integration error, as detailed in Section 2.4.1 on page 28. This leads to a maximal error in O(log(n) s-1 /n), when using a QMC estimator with n samples. In this section, we present a list of existing samplers whose goal is to ensure the best possible uniformity.

The issue with most those methods is that the resulting sampling patterns present a very regular structure and are therefore very prone to aliasing. Furthermore, as they are deterministic, it becomes difficult to decorrelate the samples between pixels. A few methods exist to address both those issues. Those algorithms are referred to as scrambling.

Here, we first present existing lattices samplers, before presenting common digital nets and digital sequences. We will finish by detailing some of the scrambling methods.

Lattices

Lattices are usually avoided in Monte Carlo rendering as there is no sampling pattern that generates as much aliasing as regular grids. However, such structures present no energy at all outside the spectral peaks, therefore, for particular integrands with limited frequencies, they can be useful. A lattice point set is defined as Definition 1 ([Lem09], Def 5.2) For a given dimension s, a lattice point set P n is defined by an integration lattice L s of the form

L s := {v 1 w 1 + ... + v s w s , v ∈ Z s } (3.12)
where the s vectors w 1 , • • • , w s in R s -which form a basis -are linearly independent over the rational numbers and are such that Z s ⊂ L s . The corresponding point set is obtained as

P n := L s ∩ [0, 1) s . (3.13)
In other words, the points in P n are obtained by taking all integer linear combinations of the vectors that fall in [0, 1) s . The resulting number of points n in P n can be shown to be equal to 1/| det(W )|, where W is the s × s matrix whose i th row is w i . To reduce the number of possible bases, we can rely on the rank r of a lattice. A Rank-1 lattice is defined as

P n := i n (z 1 , • • • , z s ) mod 1, i = 0, • • • , n -1 (3.14)
where z i ∈ N. The vector z is referred to as the generating vector. An interesting property of such lattices is that they can be made fully projection regular simply by choosing the values for z i to be prime with n.

One example of construction for the vector z is Korobov point sets [START_REF] Nm Korobov | The approximate computation of multiple integrals[END_REF] (Fig. 3.29).

A Korobov point set P n is created from a generator value a, and the following z function,

z := (1, a, a 2 mod n, • • • , a s-1 mod n) (3.15) P n := i n z mod 1, i = 0, • • • , n -1 . (3.16)
When a and n are prime relatively to each other, Korobov point sets are fully projection regular.

Dammertz and Keller introduced those lattices for rendering purposes [START_REF] Dammertz | Image synthesis by rank-1 lattices[END_REF]. Their idea is to perform an exhaustive search over the Korobov generators to maximize the minimum distance between samples of the grid. This increases the size of the zero region in the Fourier Spectrum.

As such Rank 1 lattices are very efficient to compute, such sampling pattern can prove very useful if we can ensure that there are no frequencies above this region.

For all previously presented grids, the computation of the samples requires a prior knowledge of the total number of samples to generate. To overcome this limitation, Hickernell introduces extensible lattice sequences [START_REF] Fred | Computing multivariate normal probabilities using rank-1 lattice sequences[END_REF]. The i th sample of such sequence is generated from a vector z with the following function

v i := φ b (i -1)z mod 1, i ≥ 1 (3.17)
where φ b is the radical inverse function in base b. In most cases, the z vector is found through exhaustive computer search [Hic+00; GL07].

However, even with an optimal z vector, we can point that in a typical rendered scene, we observe a very wide range of frequencies. Therefore, in most situations, the grid structure of the above mentioned methods will create massive aliasing artefacts. Instead, we will rather use other low discrepancy sampling patterns, the digital nets and sequences.

Digital Nets and Sequences

Before diving into the details of the algorithms generating digital nets and digital sequences, we will see more precisely what exactly is a digital net, and what is a digital sequence. The name digital net refers to all s-D point sets in the domain [0, 1) s , defined in a basis b, that are (t, k, s)-nets (see Definitions 2 and 3). Roughly speaking, such point sets contain b k samples, and if you partition their domain into intervals that are fractions of b q , you will have the same amount of samples inside each interval. When q = 0, there is a single interval containing all samples. The smaller the factor t = kq, (with q the highest possible value inducing a valid partitioning) the better the uniformity of the net. More formally, Lemieux [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF] derives the definitions of (t, k, s)-nets from the concept of (q 1 , • • • , q s )-equidistribution.

Definition 2 ([Lem09, Def 5.6]) Let q 1 , • • • , q s be non-negative integers, and let

q 1 + • • •+ q s . A point set P n with n = b k points is (q 1 , • • • , q s )-equidistributed in base b if every cell (or elementary interval) of the form J (r) := s j=1 r j b qj , r j + 1 b qj , (3.18) for 0 < r j < b qj , j = 1, • • • , s, contains b k-q points from P n .
The intervals J (r) are called dyadic intervals.

Definition 3 ([Lem09, Def 5.7])

A set P n containing n = b k points is called a (t, k, s)-net in base b if it is (q 1 , • • • , q s )-equidistributed in base b whenever q ≤ k -t.
Figure 3.30 on the following page illustrates a (0, 4, 2)-net in base 2. This set contains 2 4 = 16 samples, that are equidistributed for all intervals which sizes are in 1 2 q 0 × 1 2 q 1 with q 0 + q 1 ≤ 4. Figure 3.31 on page 61 illustrates a (2, 4, 2)-net in base 2. This set contains 2 4 = 16 samples, that are equidistributed for all intervals which sizes are in 1 2 q 0 × 1 2 q 1 with q 0 + q 1 ≤ 4 -2. Note that, as stated in [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF], (0, m, s)-nets in base b = 2 can only exist for s ≤ 3. This leads us to the concept of digital sequence. A digital sequence refers to point sets of infinite size. This means that from a digital sequence, we can generate a point set P n with n samples, and then, with the same algorithm, enrich this point set with m new samples to generate a new set P n+m . Usually, such sequence is defined implicitly. A digital sequence is a (t, s)-sequence, which is an extension of a (t, k, s)-net.

Low discrepancy samplers

q 0 + q 1 = 4, 2 4-4 = 1 sample per box q 1 = 0, q 2 = 4 q 1 = 1, q 2 = 3 q 1 = 2, q 2 = 2 q 1 = 3, q 2 = 1 q 1 = 4, q 2 = 0 q 0 + q 1 = 3, 2 4-3 = 2 samples per box q 1 = 3, q 2 = 0 q 1 = 2, q 2 = 1 q 1 = 1, q 2 = 2 q 1 = 0, q 2 = 3 q 0 + q 1 = 2, 2 4-2 = 4 samples per box q 1 = 2, q 2 = 0 q 1 = 1, q 2 = 1 q 1 = 0, q 2 = 2 q 0 + q 1 = 1, 2 4-1 = 8 samples per box q 1 = 2, q 2 = 0 q 1 = 0, q 2 = 1 q 0 + q 1 = 0, 2 4-0 = 16 samples per box q 1 = 0, q 2 = 0 q 0 + q 1 = 4, 2 4-4 = 1 sample per box q 1 = 0, q 2 = 4 q 1 = 1, q 2 = 3 q 1 = 2, q 2 = 2 q 1 = 3, q 2 = 1 q 1 = 4, q 2 = 0 q 0 + q 1 = 3, 2 4-3 = 2 samples per box

q 1 = 3, q 2 = 0 q 1 = 2, q 2 = 1 q 1 = 1, q 2 = 2 q 1 = 0, q 2 = 3 q 0 + q 1 = 2, 2 4-2 = 4 samples per box q 1 =
2, q 2 = 0 q 1 = 1, q 2 = 1 q 1 = 0, q 2 = 2 q 0 + q 1 = 1, 2 4-1 = 8 samples per box q 1 = 2, q 2 = 0 q 1 = 0, q 2 = 1 q 0 + q 1 = 0, 2 4-0 = 16 samples per box q 1 = 0, q 2 = 0

(l+1)b k-1 , • • • , v (l+1)b k -1 with k > t and l > 0 is a (t, k, s)-net in base b.
Furthermore, if a sequence is a (t, s)-sequence, we have the following upper bound on its l ∞ star discrepancy ([Nie88, Eq 3.])

D * ≤ c s (log(n)) s + O log(n) s-1 n (3.19)
where c s is a constant depending on the basis b, the number of dimensions s, and most important, on the t value. Formally, we have ([Lem09, Eq 5.6]) For the (t, ks)-net definition, the t value is the smallest

c s := 1 s! b t b -1 2 b/2 b/2 log(b) s . ( 3
t = q 0 + • • • + q s-1 for which P n is q 0 , • • • , q s-1 -equidistributed.
Similarly, for a (t, s) sequence, the t value is the smallest such that the subsets of P n are (t, k, s)-nets.

Now that the concepts of digital nets and sequences are defined, we will see the existing methods to generate such samplers. What definition 4 means, is that there exist an implicit link between the discrepancy and the net properties of a sampler. Even though the concept of low discrepancy is defined asymptotically and thus cannot be used to describe a single point set, we can however say that any set that is a (t, k, s)-net has a guaranteed uniformity, as good as its t value is low. Note that as digital sequences guarantee net properties for a infinite number of samples, we often refer to them as Low discrepancy sequence (LDS).

Many digital nets or sequences follow the same construction. The coordinate in each dimension d of each sample can be generated from an analytic function f d , taking as parameter the index of the sample in the set. A sample x, the i th sample of a s-D set, can be represented as

x (i) := (f 1 (i), ..., f s (i)). (3.21)
The function f d (i) calculates the coordinate in the d th dimension of the i th sample of the set. Thus, the values of the f d functions are values in the domain of the sampling pattern which is usually [0, 1) s . In most situations, the f d function computes the coordinate of the i th sample by performing bitwise modifications to the index i.

As was detailed Section 1.4 on page 8, the bit representation in base b of a value x d ∈ [0, 1) s is as follows:

x d = ∞ i=1 a i b -i (3.22)
also represented as

x d = (.a 1 a 2 • • • a N ) b (3.23)
where a k is the k th digit of the b-ary representation of x d .

Following this, the function f d can be noted as

f d (n) := C d ⊗ (a 0 , a 1 , • • • , a N ) T (3.24)
where C d is a N × N matrix, and a k is the k th digit of the b-ary representation of the index n, with N the total number of digits needed to represent n in base b. Note that as the number of samples increases, the value of n increases as well, and therefore the value of N increases too as more digits will be needed to represent a higher n in base b. Note that if we want our resulting sampling pattern to be a (0, k, s)-net, the matrices C d must follow some constraints:

Theorem 1 ([Nie92, Theorem 4.28]) Let P n be a point set of size n = 2 k in base 2 and dimension s generated from the matrices C 1 , • • • , C s . P n is a (0, k, s) net in base 2 if for all d = (d 1 , • • • , d n ) T ∈ N s 0 with |d| = m the following holds det(C s j=1 dj ) = 0 (3.25)
where

C s j=1 dj := ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ c 1 1,1 • • • c 1 1,m . . . c 1 d1,1 • • • c 1 d1,m . . . c s 1,1 • • • c s 1,m . . . c s ds,1 • • • c s ds,m ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ∈ F m×m 2 (3.26) with C j := c j k,l m k,l=1 for j = 1, • • • , s. This means C s j=1 dj is an m × m matrix consisting of the first d j rows of the generator C j for all j = 1, • • • , s.
A careful reader will have noted that Theorem 1 only applies to digital nets. Further constraints apply for generating digital sequences. We would like to first point out that if a matrix C d is not upper triangular, the resulting sampler is not a (0, s)-sequence. To be a sequence, a sampler has to be able to increase the size of the C d matrices on the fly, without changing the result of the computation of the previous indices. Therefore, all values below the diagonal must be 0, to act as neutral elements in the computations and therefore preserve the position of the previous samples. When generating digital nets, we do not need to increase the size of the matrix on the fly as we already know the number of samples required. Therefore, we also do not need the matrix to be upper triangular. We will detail this further as we present the common existing digital sequences and digital nets.

Digital Sequences

The first and simplest digital sequence is the Van der Corput sequence [Cor35]. This 1-D sequence generates the coordinates of the i th sample by reversing its digit in base 2. This translates in Equation 3.24 on the previous page as This sequence was extended by Halton [START_REF] John H Halton | Algorithm 247: Radical-inverse quasi-random point sequence[END_REF] to a s-D sequence through the use of different bases in each dimension. When you perform a 1-D bit reversal in a basis b, it sums up as dividing your domain into b n intervals and putting one sample inside each interval. In s-D, if the bases used are prime numbers on all dimensions, one is guaranteed that the intervals will not be aligned between dimensions and therefore the sequence will be low discrepancy. Furthermore, this use of several basis also helps reducing the regularities and therefore the peaks in the resulting Fourier spectrum (Fig. 3.32).

(a 0 , a 1 , • • • , a N ) T = I ⊗ (b 0 , b 1 , • • • , b N ) T (3.

The Sobol Sequence

For sequences using a single basis, the most famous (t, s)-sequences is probably the Sobol sequence [Sob67] (Fig. 3.33). Defined in base 2, this sequence builds the matrices C d from primitive polynomials. The choice of the primitive polynomials used will affect greatly the result of the sampling pattern (Fig. 3.34 on the following page).

The list of all the primitive polynomials usable for this sequence can be found in [Oei]. Each polynomial is referred to using its index in this list. We build each C d matrix from the polynomials using the following algorithm, which can also be found in [START_REF] Joe | Notes on generating Sobol sequences[END_REF]. We remind that each dimension uses a different C d matrix, generated from its corresponding primitive polynomial.

We start with a primitive polynomials of degree k in the field F 2 , denoted as

x k ⊕ a 1 x k-1 ⊕ a 2 x k-2 ⊕ • • • ⊕ a k-1 x + 1 (3.28)
where the coefficients a i are either 0 or 1. A set of numbers m 1 , m 2 , • • • is built from the following recurrence relation

m j := 2a 1 m j-1 ⊕ 2 2 a 2 m j-2 ⊕ • • • ⊕ 2 k-1 a k-1 m j-k+1 ⊕ 2 k m j-k ⊕ m j-k (3.29)
You may have noticed from this equation that it requires an initial subset of m j value to iterate upon. Those values can be chosen freely, as long as for each m j with 1 ≤ j ≤ k, m j is odd and m j is lesser than 2 j . Those conditions ensure that the binary C d matrix formed for the m j values will present 1 on its diagonal and will be upper triangular. Those numbers are used to create the set of direction numbers Finally, the C d matrix is created by concatenating the v j values, expressed in binary, as its columns. Meaning, if we denote v i j the i th bit of the number v j , we have

{v 1 , v 2 , • • • } with v j = m j 2 j (3.30)
C d = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ v 0 1 v 0 2 v 0 3 • • • v 1 1 v 1 2 v 1 3 • • • v 2 1 v 2 2 v 2 3 • • • . . . . . . . . . . . . ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 v 0 2 v 0 3 • • • 0 1 v 1 3 • • • 0 0 1 • • • . . . . . . . . . . . . ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
(3.31)

The simplification to an upper triangular matrix is guaranteed from the conditions on the initial m j values.

The Sobol sequence has been proven to be a (t, s)-sequence. However, the higher the degree of the polynomials, the higher the t value and therefore the worse the quality of the sampling pattern. Formally, we have

t := s d=1 k d -1, (3.32)
where k d is the degree of the polynomial used in the d th dimension [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF]. It is interesting to know that there is a limited number of primitive polynomials of each degree, and that to be a (t, s)-sequence, the Sobol sequence can not use the same polynomial twice. Therefore, the higher the dimension of the point set, the worse its uniformity. This is illustrated in Figure 3.34.

This drawback is actually true for most of the low discrepancy sequences. One way to alleviate this issue is to permute the digits of the sample index before computing its coordinate. This does not change the position of the samples of the sequence but it changes their order in this sequence. This can be done using Gray Code. Each index i is replaced with its Gray Code equivalent, denoted G, as

G(i) := i ⊕ i 2 (3.33)
where the ⊕ operator is the logical xor operator, and [.] represents the integral part function.

Note that the function G takes as input an integer i ∈ N, and outputs another integer G(i) ∈ N.

P 16
The first 12 original samples

The first 12 Gray code samples Applying Gray Code changes Equation 3.24 on page 63, that generates the d th dimension of the i th sample of the s-D sequence defined from a set of s generative matrices

{C d }, to (c N , • • • , c 0 ) 2 = G(i)(a 0 , a 1 , • • • , a N ) = C d • (c 0 , c 1 , • • • , c N ) T . (3.34)
The impact of Gray Code is illustrated Figure 3.35.

The Sobol sequence was later extended in several ways. In [START_REF] Tezuka | Uniform random numbers: Theory and practice[END_REF], Tezuka obtained a sequence no longer restricted to primitive polynomials. In [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF], the Faure sequence was developed. This sequence aims at building (0, s)-sequences in base b, relying on Pascal Matrices raised to the power j, j = 0, • • • , s. Note that this construction requires that s ≤ b, which prevent a dynamic increase of the number of dimensions. Finally, Niederreiter [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF], extended this sequence to bases b = 2, which were further extended to Niederreiter-Xing sequences [START_REF] Niederreiter | Low-discrepancy sequences obtained from algebraic function fields over finite fields[END_REF], where t increases linearly with s (instead of exponentially for other sequences), which is the optimal achievable rate.

All the sequences presented above allow for a direct generation of the i th sample, without any prior knowledge of the total number of sample of the point set. This property adds constraints on the sample generation that limit the number of possible sampling patterns. If we choose to remove the sequentiality, we will generate digital nets, which are less constrained than digital sequences and therefore allow for new sampling patterns.

Digital Nets

An interesting property of digital sequences is that any digital sequence defined in dimension s can be turned into a digital net in dimension (s + 1) [START_REF] Niederreiter | Random Number Generation and quasi-Monte Carlo Methods[END_REF]. To do so, simply set the matrix C s+1 as

C s+1 := ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 • • • 0 1 0 0 • • • 1 0 0 1 • • • 0 0 1 0 • • • 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (3.35)
This can be intuitively defined as setting the (s + 1) th coordinate of the i th sample in a point set of size n as i n .

The Hammersley point set [START_REF] John | Monte Carlo methods for solving multivariable problems[END_REF] is the combination of this property with the Van der Corput sequence (Fig. 3.36 on the facing page). This set is defined in 2-D, where the i th sample x is defined as

x := φ(i), i n (3.36)
with φ being the Van der corput sequence and n being the size of the point set. Expressed in the framework of Equation 3.24 on page 63, we have

C 1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 • • • 0 0 0 1 • • • 0 0 . . . 0 0 • • • 1 0 0 0 • • • 0 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (3.37) C 2 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 • • • 0 1 0 0 • • • 1 0 . . . 0 1 • • • 0 0 1 0 • • • 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (3.38)
Due to the high structure induced by the implicit definition of the sampling patterns, Digital nets are deterministic and regular sampling pattern. When used in rendering, this generates aliasing and correlation between pixels. To alleviate this, scrambling algorithms were devised and will be presented in Chapter 4.

Conclusion

In conclusion, among existing samples we can denote three major families. The stochastic historical samplers, the Blue Noise samplers, and the Low Discrepancy analytically defined samplers. Each of those families have its pros and cons, as illustrated on the graphs 3.37 and 3.38.

On the graph 3.37, we roughly compare the performances of the main samplers for each family, with respect to their spectrum quality, and their discrepancy. A spectrum is considered of good quality if it is close to the Blue Noise spectrum and does not present spectral peaks. For the discrepancy, the right most samplers are those that are Low Discrepancy. We can note that the best Low discrepancy spectrum is also the one with the worst Fourier spectrum, and that the discrepancy of the sampler with the best spectrum is only average.

The graph 3.38 regards mostly the practicality of various samplers. A sampler can present the best spectrum and Low discrepancy, if it can't be used in practice it is rather pointless. On this graph, we compare samplers with regard to their computational efficiency (the speed at which they can generate high amount of samples), and their dimensionality (a higher dimensionality means a sampler can generate samplers in higher dimensions). Both those characteristics are mandatory for rendering, as we solve recursive higher-dimensional integrands, and generate several hundreds of samples per pixels. It can be seen on this graph that canonical and Low Discrepancy samplers are very practical, when Blue Noise samplers are almost impossible to use in a rendering context, despite the quality of their spectrum.

Following this, we conclude that historical samplers have a spectrum of reasonable quality (in the sense that is does not generate aliasing) and are efficient for higher dimensions. Blue Noise samplers present an optimal spectrum and are therefore ideal for Monte Carlo Rendering as they have a very low variance and do not generate aliasing. However, such samplers are limited to low number of dimensions, and are not computationally efficient. Finally, the low discrepancy samplers present a very low error in integration, are very efficient even in high dimensions thanks to their analytical definition, and some also allow for iterative addition of samples. However, their structured sampling patterns induce strong aliasing when used in Monte Carlo rendering. Table 3.1 on the following page sums all the integration errors (or MSE) of various samplers, presented in this chapter.

Sampler Worst Case MSE/Error 2 Discrepancy Whitenoise O(1/n) O(1/ (n)) Stratified O(1/(n √ n)) Poisson Disk O(1/n) Blue Noise O(1/(n √ n)) Low Discrepancy O((log(n)/n) 2 ) O(log(n)/n)
Tab. 3.1.: Worst case convergence speeds of the integration error when using warious sampling patterns, along with their discrepancy.

In this PhD, we aimed at creating hybrid samplers. Those samplers would have both a Blue Noise spectrum and be Low Discrepancy. This allows them to have the lowest error in integration, without any aliasing, while still being computationally efficient in higher dimensions.

Optimizing a set

As shown in Section 2.1.1 on page 17, the discrepancy of a set creates an upper bound on the error of the Monte Carlo integration using this set. Therefore, the lower the discrepancy the lower the bound on the error. However, most low discrepancy samplers are highly structured, which generates aliasing when used for rendering an image. On the other hand, among stochastic samplers, the best distribution of samples for MSE reduction is the Blue Noise distribution (Section 2.4.2 on page 30). This distribution also does not present any spectral peaks and thus does not generate aliasing. We thus try to create sets that are both Low discrepancy and Blue Noise, in hope to reach an optimal low in integration error along with an anti aliased spectrum.

In this Chapter, we present a few naive methods to improve the spectrum or the discrepancy of various samplers. The first category of methods, Section 4.1, aimed originally at randomizing Low Discrepancy analytic samplers. However, we will see that one of those methods, named Owen's scrambling [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF], also allows to improve the Fourier spectrum of point sets.

In a second time, Section 4.2 on page 79 presents our adaptation of Owen's scrambling that allows to generate samples sequentially. Finally, Section 4.3 on page 80 and Section 4.4 on page 84 presents our early attempts to either improve the discrepancy of a Blue Noise set or to scramble a Low Discrepancy set to give it a Blue Noise Fourier spectrum.

Randomizing Low discrepancy sets

In general, all the Low Discrepancy point sets and sequences presented in Chapter 3 show very regular point placement. This leads to peaks in their Fourier spectrum which in turn lead to aliasing. Furthermore, those point sets are deterministic, which is an issue in rendering as using the same point set in each pixel also leads to aliasing due to pixel correlation. In this section, we present some common randomizing algorithms to alleviate those correlations.

We first present trivial shifting algorithms, which do not impect the spectral properties of the final set. Thus, we will the nmove on to presenting a specific 1-D scrambling method, that we will reuse in Chapter 5 on page 93 to adapt the spectrum of a Low Discrepancy set. Finally, we will present Owen's scrambling, a n-D scrambling that will be reused in Chapter 6 on page 107.

Shifting

The easiest method to randomize point sets are the digital shifts. Such methods combine each sample of the set with a single randomly chosen input. They thus do not break the sequentiality of an input sequence.

The most straightforward one is the Cranley-Patterson shift [START_REF] Cranley | Randomization of number theoretic methods for multiple integration[END_REF]. This applies a toric digital shift, chosen randomly, to a point set. More formally, with this method, a sample x is shifted to create a new sampler x from a shift vector v with the following formula

x := x + v mod 1, (4.1)
where the + operator performs the coordinate wise sum of both vectors.

This scrambling efficiently decorrelates the pixels in rendering while preserving the Fourier spectrum of the set. Following [KN74, Lemma 1.1], such method preserves the uniformity of the given set. Note that if this scrambling is applied on a Digital net, it does not preserve the net properties of the input set.

If one needs to preserve the net properties of the set, one can use the Digital shift. This method takes a vector at random and performs a xor in base b between this vector and all samples. If we denote x the Digital shifting of x with a vector v, we have

x := x ⊕ v, (4.2)
where the ⊕ operator performs the coordinate wise xor of the two vectors.

Those trivial shifting methods are extremely efficient when it comes to randomizing a point set. However, they pretty much do not do anything else. The following sections will therefore present other randomizing methods, that affect the Fourier Spectrum of a point set.

Chunk based 1-D scrambling

In [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF], we can find the following theorem,

Theorem 2 ([KN74, Chapter 2, Theorem 4.1]) Let P n := {x 1 , x 2 , • • • x n } and P n := {y 1 , y 2 , • • • , y n } be two finite sequences in [0, 1). Suppose 1 , 2 , • • • , n are nonnegative numbers such that |x k -y k | ≤ k for 1 ≤ k ≤ n.
Then, for any ≥ 0, we have

|D * (P n ) -D * (P n )| ≤ 2 + N ( ) n (4.3)
where N ( ) denotes the number of k, 1 ≤ k ≤ N , such that k > .

Randomizing Low discrepancy sets

This theorem starts with two 1-D indexed sets of samples, such that for all i, the i th sample of the first set can be paired with the i th sample from the other set. The theorem then states that, if the difference between paired samples does not exceed a value , then the difference between the discrepancies of both input sets is in O(1/n), and is bounded by 2 + n/n.

In our case, we will use this theorem to evaluate how much scrambling a 1-D Low discrepancy sequence affects the discrepancy of the original sampling pattern. Instead of having two different sequences as input, we will take a 1-D Low Discrepancy sequence S, and then its scrambled version. We then see that if we scramble it locally, the difference between the discrepancy of those sequences is bounded by the size of the scrambled subset. We will come back to this in Chapter 5 on page 93, where we will give formal definitions and formal proofs for this. We will also see how this 1-D theorem can be used in a 2-D context for spectral control of Low Discrepancy sequences.

Owen's scrambling

There exists another method, that scrambles sets of higher dimensions, Owen's scrambling [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF], (already used for Rendering in [START_REF] Veach | Robust monte carlo methods for light transport simulation[END_REF]). This method allows to preserve the discrepancy of a full s-D set by randomly scrambling each 1-D projections of this set. It can also impact the Fourier Spectrum of the initial set, and therefore served as a basis for most our experiments.

For the particular case of digital nets, advanced scrambling method exists [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF]. The most important one for us is the Owen's scrambling [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF], also referred to as fully random scrambling.

This scrambling usually applies on sets that are (t, k, s)-nets, and guarantees to preserve the net properties of the set. It could be applied on other categories of sampling patterns but would not be as easy to implement and we would not have any guarantees on the resulting set.

It relies on swapping some coordinates of some samples, randomly chosen. This swapping is performed independently over each dimension. A s-D point set is permuted using s 1-D permutations, one for each dimension. Therefore, we can detail this scrambling in 1-D without any loss of generality.

Geometrically, a domain Ω is recursively subdivided by a factor 2. Then, a binary tree of boolean flags is used to determine if two subdivisions of space are to be swapped or not. Each flag of this tree corresponds to a pair of neighbouring subsets of the domain, and by filling this tree of boolean flag with random values, we will randomly swap portions of space. This is illustrated in 1-D Figure 4.1 on the next page and in 2-D Figure 4.2 on the facing page.

We now give a more formal definition of this scrambling. For a set of n 1-D samples P n , if P n is a (t, k, 1)-net, each sample can be expressed as an integer, defined in base 2 by d = log 2 (n) bits. Thus, Owen's scrambling of this set requires a binary tree of depth d. This tree is filled with a set of n -1 boolean flags chosen randomly. Formally, a boolean tree is defined as a set of flags {B j }, where j is the binary vector representing the id of this flag at a given depth i.

The length of the binary vector j allows to retrieve the level i of a given flag. This boolean tree is then applied on each 1-D sample of the set. Let us consider a 1-D sample with coordinate x ∈ [0, 1) defined in base 2 with digits a i (i.e. x := ∞ i=1 a i 2 -i ). Owen's scrambling creates a new 1-D point π(x) from a tree of depth d. The initial sample x is defined in base 2 as x = (.a 0 a 1 • • • a N ) Each digit a k of this point is xored with the flags in the permutation tree as follows, leading to a new set of b k digits:

B (.) B (0) B (
b 0 := B (.) ⊕ a 0 b 1 := B (a0) ⊕ a 1 b 2 := B (a0a1) ⊕ a 2 . . . b d-1 := B (a0a1•••a d-2 ) ⊕ a d-1
where the ⊕ operator performs the binary xor between two boolean values. This leads to a final scrambled point

π(x) := (.b 0 b 1 • • • b d-1 a d a d+1 • • • a N ) (4.4) (note that if N < d, we have ∀k, k > N, a k = 0).
This means that Owen's scrambling goes through all digits and swaps the digit a i if its associated flag is set to 1 in the tree at depth i. Geometrically, swapping a digit a i swaps all samples with the same a 0 ...a i binary prefix (thus all the samples that geometrically belong to the same portion of space). Owen's scrambling in 1-D therefore implies successive permutations between intervals of points of size 2 -i+1 . Please refer to [START_REF] Grünschloß | Enumerating quasi-monte carlo point sequences in elementary intervals[END_REF] for efficient implementations of this algorithm.

Note that this scrambling, when applied on a digital net, has the interesting property of both preserving the Low Discrepancy property of this net (meaning that it preserves its equidistribution), and of being able in some cases to remove the peaks in its Fourier Spectrum 4.3. This removal of the peaks happens with a probability one when the permutation trees, that control if two subintervals of space are swapped, are filled at random with a high enough depth d. This leads to a random exchange of samples coordinates, which breaks the original structure or the net.

4.2 Owen's scrambling, from a random set to a random sequence

In the previous section, we presented the original Owen's scrambling algorithm. This algorithm manages sets of samples. As the depth d of the tree is fixed, this algorithm implicitly assumes that we have only 2 d samples, and therefore, it is not sequential. If one tries to iteratively increase the size of this tree with randomly chosen flags, maintaining the consistency between samples would require re-evaluating this tree for all samples, which might move previously scrambled samples, thus breaking the sequentiality.

Our goal is to be able to iteratively increase this tree without having to re-evaluate this tree for previously existing samples, while still preserving the consistency of the whole sequence. To do so, we rely on the two following observations:

• Each boolean tree of depth d has 2 d-1 different branches, identified as the binary vector of the id of their leaf node (a 0 , a 1 , . . . , a d-1 ). Each branch applies to the pair of samples whose d -1 first digits are (a 0 , a 1 , . . . , a d-2 ).

• If all the flags along the branch (a 0 , a 1 , . . . , a d-1 ) are 0 values, the samples with coordinate (.a 0 a 1 . . . a d-2 0) 2 and (.a 0 a 1 . . . a d-2 1) 2 will not be moved by scrambling.

Following those observations, it appears that we can iteratively increase the size of a permutation tree, creating a new tree of depth d + m, without affecting previously existing samples. To do 4. [START_REF]D point set with stratified 1-D projections[END_REF] Owen's scrambling, from a random set to a random sequence so, we need to figure out which branches apply to those samples, and to fill the extensions of those branches with flags at 0.

For example, the following tree of depth 3:

B (.) B (0) B (1)
B (00) B (01) B (10) B (11)
can be extended, while preserving sequentiality, to the following tree of depth 5:

B (.) B (0) B (1) B (00) B (01) B (10) B (11) 0 B (••• ) 0 B (••• ) 0 B (••• ) 0 B (••• ) 0 B (••• ) B (••• ) B (••• ) 0 B (••• ) B (••• ) B (••• ) 0 B (••• ) B (••• ) B (••• ) 0 B (••• ) B (••• ) B (••• )
Depending on the size of m, only m • 2 d flags are required to be 0, leaving the remaining 2 d+mm • 2 d flags to be randomly chosen. Note the trade-off regarding the m value, if m is too small, the point set will not be properly anti-aliased, (Fig. 4.4 on the next page) if m is too big, you may end up will evaluating a very deep tree for only a small set of samples.

With this contribution, one can use Owen's scrambling in a high dimensional sequential context. However, it still does not generate sets that are Blue Noise. A further extension of Owen's scrambling is presented in Section 4.4 on page 84 that allows generating small patterns with Blue Noise spectrum and Low discrepancy. Before that, we present in Section 4.3 another contribution that aims at improving the discrepancy of a stochastic pattern.

From an Anti-aliased to a Low Discrepancy set

In this section, we present a simple method that can be used to reduce the discrepancy of a pre-existing set. In its current state, our method takes as input a given set (for example from the BNOT sampler [START_REF] De Goes | Blue noise through optimal transport[END_REF]), and optimizes its generalized l 2 star discrepancy. Where the l ∞ star discrepancy keeps the maximal observed difference, the generalized l 2 star discrepancy averages all the squared differences observed from the area of anchored boxes and samples:

(D 2 (P n )) 2 := 4 3 s - 2 n n i=1 s j=1 3 -x (i) j 2 2 + 1 n 2 n i=1 n i =1 s j=1 2 -max x (i) j , x (i ) j . (4.5)
The main idea, in this section, is to solve the following optimization problem:

Find point set P * (n) such that D 2 (P * n ) is minimal.
We follow a variational approach expressing the variation of the discrepancy as we slightly move a given point (i.e. its derivative with respect to sample positions). First, we observe that max terms in Eq. (4.5) cannot be differentiated. We first rewrite (4.5) using a smooth approximation of the max (sometimes referred as softmax):

M α (u, v) := ue αu + ve αv e αu + e αv . ( 4.6)
This smooth approximation tends to the maximum value between u and v as α increases. A key property is that M can be differentiated with respect to u (similarly to v):

∂M α ∂u (u, v) = e α(u+v) (α(u -v) + 1) + e 2αu (e αu + e αv ) 2 . ( 4.7)
Replacing the max by M α in (4.5) leads to a differentiable definition of the (square of the) generalized l 2 discrepancy. More precisely, for a given 2-D sample x (i) := (x i , y i ), we obtain the explicit gradients of discrepancy at this given point along its components:

∂(D 2 (P n )) 2 ∂x i = -2 + y i + 2 n j=1
2 -e αx i + e αx i αx i e αx i + e αx j + αe αx i (e αx i x i e αx j x j ) (e αx i + e αx j ) 2 2 -e αy i y i + e αy j y j e αy i + e αy j , (4.8)

∂(D 2 (P n )) 2 ∂y i = -2 + x i + 2 n j=1
2 -e αy i + e αy i αy i e αy i + e αy j + αe αy i (e αy i y i e αy j y j ) (e αy i + e αy j ) 2 2 -e αx i x i + e αx j x j e αx i + e αx j .

(4.9)

We can now combine all independent gradients into a single vector

∇(D 2 (P n )) 2 := ∂(D 2 (P n )) 2 ∂x i , ∂(D 2 (P n )) 2 ∂y i i=0...n-1
defined for all points i of P n . For convex energy function, optimal sample positions P * n would have been obtained solving ∇(D 2 (P n )) 2 = 0. In our context, Eq.(4.5) is not convex and has many local minima. We minimize the non linear energy using a simple Gradient Descent approach as described in Algorithm 1: we move each point of a given pattern to the opposite direction of its gradient by a factor (until the gradient norm is below a given threshold). Many improvements exist using for instance a line search to find a better λ at each step, or higher order schemes using Newton's or Quasi-Newton's approaches. We just describe here a proof of concept that already gives interesting results. Again, the fixed point of the algorithm may not be a global minimum but the improvement in terms of discrepancy is valuable while preserving the spectral properties of the starting point set.

Algorithm 1: Simple gradient descent algorithm to minimize the generalized l 2 star discrepancy.

Input: A 2-D point set P n to optimize and a threshold τ .

Output: A 2-D point set P n with a lower discrepancy.

1 repeat 2 l 2 norm ← 0; 3 forall x (i) ∈ P n do 4 ∇x (i) ← computeGradient(x (i) ); 5 x (i) ← x (i) -λ • ∇x (i) ; 6 l 2 norm ← l 2 norm + ∇x (i) 2 ; 7 end 8 until l 2 norm < τ ; 9 return P n ;
To compute the gradient at a point x (i) (line 4), we have experimented several strategies: Explicit smooth gradients: use explicit gradients using Eq. (4.8) and (4.9). As discussed earlier, we need the smoothing factor α to be large to have a robust approximation of the max operator. The main issue is that as the number of points increases, we need to also increase α to refine the approximation. For large point sets (n > 500), we observe numerical inconsistencies leading to incorrect results or non-convergence of Algorithm 1.

Auto-differentiation of the non-smooth generalized discrepancy: another solution consists in considering the non-smooth original definition of the Generalized discrepancy (hence with the max). In that case, the gradient is obtained using automatic differentiation of the function computing the energy. The auto-differentiation tool ( [Dcp]) will differentiate the singularities around the "if" statements involved in the energy expressions for the max terms, leading to an expression of the code-based gradient that we can use in the Gradient Descent1 .

Results given below have been obtained using the second approach. , and the Fourier spectrum of this same set after optimization. We can note that a black cross appears in the middle of the spectrum. This cross appears as, to reduce the discrepancy, the samples aligned on either axis were pushed away from each other by our optimization. For now, this method is still limited but we hope that, by combining it with Oztirelli's method to fit a point set to a given PCF [START_REF] Öztireli | Analysis and Synthesis of Point Distributions based on Pair Correlation[END_REF], and by porting on the GPU, we could generate very quickly stochastic low discrepancy sets with a given Fourier spectrum, even for a high number of samples. We further could add an adaptivity criterion to our sets, using local discrepancies measures and PCF. But this is all for future work.

From a Low Discrepancy to an Anti-aliased set

The methods presented in this section aim at scrambling/generating Low discrepancy sets to control their Fourier spectrum. We start by presenting an existing method, that optimizes generative matrices from exhaustive search (Section 4.4.1 on page 87). We then present our new contributions, based on genetic algorithms, to generate Blue Noise sets from Owen's scrambling on the Sobol sequence (Section 4.4.2 on page 87).

Generative matrices

To generate Low discrepancy sets, a solution is to generate digital nets in base b. To do so, one needs s b-ary matrices (one for each dimensions), constrained by the Theorem 1 on page 63. Then, by computing the matrix vector product of the b vector of the index of samples and the appropriate matrix for each dimension, one generates a set of s-D samples that form a digital set.

In [START_REF] Grünschloß | (t, m, s)-Nets and Maximized Minimum Distance[END_REF], authors aim at generating 2-D digital sets while maximizing the minimal distance between samples. To do so, they start by exploiting a previously mentioned property of s-D digital sequences; they can be turned into (s+1)-D digital sets by setting one of their generating matrices as a mirror identity matrix. In their paper, the authors want to generate (0, k, 2)-nets maximizing the minimal distance between pairs of samples of the set. They do this by setting the matrix C 1 as a mirror of the identity matrix:

C 1 := ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 • • • 0 1 0 0 • • • 1 0 . . . 0 1 • • • 0 0 1 0 • • • 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (4.10)
Then, they simply perform an exhaustive search on all the possible generating matrices to find C 0 , and keep the matrices that generated sets with the maximal minimal distance. This kind of exhaustive search limits them to 2-D sets of small size (up to 1024 samples), therefore, they suggest that those matrices can be increased by padding them with the identity matrix.

Even though they did manage to have a higher minimal distance between samples, we cannot say that their spectrum is Blue Noise or aliasing free (Fig. 4.8 on the following page). Their issue is that, to make their exhaustive search feasible, they enforce the C 1 matrix, and, as that they generate digital nets, the matrix C 0 is highly constrained. Thus, they do not have a lot of freedom to generate their sets of samples.

We tried to alleviate some of those issues in Section 4.4.2 by using genetic algorithms and Owen's scrambling to perform our searches.

Genetic algorithms

Our goal is to optimize a Low Discrepancy point set. To do so, one way is to rely on a definition of this point set using a numerical structure, that can only take a finite number of values.

In the case of [START_REF] Grünschloß | (t, m, s)-Nets and Maximized Minimum Distance[END_REF], their point sets were defined by a binary matrix. In our contribution, we chose to define our sets as a set of Owen's permutation trees, along with the initial Sobol set to scramble (each distinct initial set requires distinct permutation trees). Note that for a s-D point set of size n, we have 2 sn possibilities for those permutation trees.

We can thus generate optimal patterns for 2-D sets with a small number n of samples through an exhaustive search on all the possible Owen permutation trees. But, for merely n = 16 and s = 2, we need to test 2 32 possibilities.

This exhaustive search is therefore impossible for higher dimensions or higher n values. Fortunately, the issue of finding the optimal solution in a space way too big to be entirely covered is well known in many domains, there are many algorithms that can optimize our search. We chose here to use genetic algorithms.

Our method initializes a population of permutation trees with random values. Each of those trees, once applied on the input point set leads to a scrambled set, that we will denote as an individual. Then, we sort those individuals using either the distance between their PCF and a target PCF, or simply by computing the minimal distance between samples (for Blue Noise).

Once those individuals are sorted, we retrieve the best ones, that will form a set denoted the elite set. We then create new individuals by either crossing our elite individuals, or generating new random values. In our experiments, we crossed individuals by exchanging chunks of their permutation trees. We then keep iterating like this a finite number of times, or until we reached a good enough individual. This algorithm is given as Algorithm 2 on the facing page.

Note that our method can become very expensive as the size of the population and the number iterations grows. This expensiveness of the method is partly alleviated as trees can be randomized. To randomized an optimized tree of depth d, simply keep the first d -1 levels as they are and fill the d th level with random values. The movement on samples will be minimal and will not affect much the quality of the final spectrum.

For small pointsets, we compare the results of this algorithm with the results from an exhaustive search in Figure 4.9 on page 90. Figure 4.10 shows how our optimization evolves with the number of iterations. Note that on both those figures, only small point sets are presented. As it is impossible for such small sets to be characterized using a PCF and a Fourier spectrum, we used the minimal distance between samples, denoted D. The higher this distance, the better the point set. The results of our genetic search for sets with higher number of samples are given in Figure (Fig. 4.11 on the facing page). This figure shows that for larger point sets, we no longer observe any improvement in the spectrum. This is due to the number of solutions that grows too much for a naive algorithm such as ours to find anything worthwhile.

Both methods, our genetic algorithms and the matrices from [START_REF] Grünschloß | (t, m, s)-Nets and Maximized Minimum Distance[END_REF] have the same drawbacks. They both present results that are an improvement (though but not as good as we would like them to be), but they are so computationally expensive that they are impossible to use in practice. In the following chapters, Chapter 5 and Chapter 6, we will present more advanced methods of the same category (trying to optimize a Low Discrepancy set to generate a Blue Noise Low Discrepancy set), that generate efficiently sets or sequences with a spectrum of higher quality.

Conclusion

All the methods presented in this section aim at solving the same issue, combining an optimal integration error (as the one of low discrepancy samplers) with the anti aliased properties of stochastic (and most of all Blue Noise) Fourier spectrum. The main ideas exposed in Section 4.4 on page 84 follow the same simple idea; finding a numerical representation for low discrepancy sets or sequences and optimizing this value, either from optimization algorithms or simple brute force searching. Those methods give very good results but only for very limited number of samples and dimensions. We definitely believe that with more advanced optimization algorithms and the use of embarrassingly parallel components such as GPUs, we could get (in the case of Owen's scrambling) a very good optimization tree that could have practical applications. However, such a method is still to be developed.

Therefore, the following Chapters present instead two very specific methods to generate first, Blue Noise Low discrepancy 2-D sets (Chapter 5 on page 93) and then, a s-D sequential sampler with Low discrepancy Blue Noise 2-D projections (Chapter 6 on page 107).

Low Discrepancy Blue Noise point sets

In Chapter 4 on page 74, we saw how one can transform a Low Discrepancy set to improve its Fourier Spectrum. However, most of previous methods required exhaustive search or complex optimization schemes and are thus difficult to use in practice. Fastest algorithms, such as Owen's scrambling do improve upon the spectrum but are still far away from the ideal Blue Noise spectrum.

In this Chapter, we present a more advanced method that efficiently generates Low Discrepancy and Blue Noise 2-D point sets [START_REF] Abdalla | Low-discrepancy blue noise sampling[END_REF] (s-D sequences will be discussed in Chapter 6 on page 107). We achieved this by scrambling an input template Low Discrepancy point set to match its spectrum to the spectrum of a target stratified point set. We use Van der Corput sequence to create a stratified low discrepancy initial set and then, through coordinate wise rearrangement, we scramble it without affecting the discrepancy of the original set. In the following, Section 5.1 presents our input Low discrepancy construction. Then, Section 5.2 on page 98 presents our scrambling method.

Initial Stratified Low Discrepancy set

Our method starts with a template low discrepancy (LD) set, that will then be optimized to acquire a desired spectral profile (BNOT [START_REF] De Goes | Blue noise through optimal transport[END_REF],

Step [START_REF] Heck | Blue noise sampling with controlled aliasing[END_REF], etc). There are many LD constructions we could have chosen from as our input LD set (Sobol, Halton, Hammersley, Rank-1, etc.). The issue is that with such sets, it is difficult to identify the neighbourhoods of the points. As the spectral properties are sensitive to local neighbourhoods [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF], we need to devise a new LD construction that addresses this issue, based on stratification. In the following, we start by formally defining a 2-D stratified point set. Then, we list the conditions to create a 2-D low discrepancy stratified point set. Finally, we present our choice as the template LD construction.

Stratification

In this section, we define formally a 2-D stratified point set. Suppose that we have a stratified point set P n , with n := k 2 samples. We denote such a set as We point out that following this definition, a point set can only contain k 2 , k ∈ N samples. Furthermore, such definition generates samples in the domain [0, k) 2 . To sample the unit domain [0, 1) 2 , one can simply scale the point set by 1 k .

P n := (X + a X,Y , Y + b X,Y ) : X, Y ∈ {0 . . . k -1}; a X,Y , b X,Y ∈ [0, 1) . ( 5

LD stratified point set

Stratified point sets are not naturally low discrepancy. Therefore, we need to develop a new construction that will let us access samples from the coordinates of their stratum while being a low discrepancy set. To do so, we generate stratified point sets P n where, for every column X, we set the offsets values from a 1-D LD sequence {a X,0 , a X,1 . . . , a X,k-1 }, and similarly for every row Y , with the 1-D sequence {b 0,Y , b 1,Y , . . . , b k-1,Y }. We can prove that this new construction is low discrepancy.

Intuitively, we can observe that the difference between the area of an integer set of strata and the fractions of samples of P n that are located into those stratas is 0. Therefore, the discrepancy of the stratified set only depends on the discrepancy of the sequences used as offset in the rows and columns.

More formally, we denote by a := {a i } a 1-D sequence of points in [0, 1), a(k) denoting the set of the first k elements of the sequence. We consider a collection of k such sequences a (X) := a (X) i , indexed by an integer X; that is, a (X) is the X th sequence in the collection, while a (X) i is the i th entry in that sequence. We then consider a similar set up for b, b, and Y (b := {b i } and b(k)) respectively, which gives us the following lemma:

Lemma 1 If n = k 2 denotes the total number of points in a stratified set P n , P n being scaled to [0, 1) s , we have

D * (P n ) ≤ 1 n max X,Y ∈{0...k-1} XD * (b (X) ) + Y D * (a (Y ) ) + 1 . (5.2)
By proving Lemma 1, we will then prove that our construction generates low discrepancy sets, as long as a and b are Low discrepancy sequences.

Proof 1 We assume that P n is a 2-D point set defined over the domain [0, 1) 2 (note that there is no loss of generality w.r.t. a point set defined over the domain R 2 ). We remind that the l ∞ star discrepancy of P n is defined as

D * (P n ) := sup v∈[0,1) 2 |λ(v) -α(P n , v)| (5.3)
where λ(v) is the Lebesgue measure of the interval [0, v) 2 , and α(P n , v) returns the fraction of points of P n that belong to the interval [0, v) 2 . We will denote d(P n , A) the 1-D discrepancy of P n over a 1-D interval A. , indexed by an integer X; that is, a (X) is the Xth sequence in the collection, while a (X) i is the i-th entry in that sequence. We consider a similar setup for b, b, and Y , respectively.

Additivity of the Lebesgue measure (s-dimensional volume) and α(P n , v) leads to the following property: Let

As P n with n = k 2 is a 2-D scaled stratified point set, composed of m low discrepancy sequences it can be denoted as

P n := X + a (X) Y k , Y + b (Y ) X k ; X, Y ∈ {0 . . . k -1} . (5.5)
In other words, the abscissa X k of a given point of the set is shifted by

a (X) Y
k , whose numerator is the Y th element of the sequence a (X) . In the following, we simplify the notations of these sequences to a and b (and their points {a i } and {b i }) in expressions whenever there is no ambiguity on X and Y values.

0 1 X k X+x k X+1 k Y k Y +1 k Y +y k 1 J XY J Xy J xy J xY
J = J XY + J Xx + J xY + J xy (5.7)
where

J XY := [0, X/k) × [0, Y/k) , J Xy := [0, X/k) × [Y/k, (Y + y)/k] , J xY := [X/n, (X + x)/k] × [0, Y/k) , J xy := [X/n, (X + x)/k] × [Y, (Y + y)/k] .
Thanks to the additive property of discrepancy, we have:

d(P n , J) ≤ d(P n , J XY ) + d(P n , J Xy ) + d(P n , J xY ) + d(P n , J xy ) . (5.8) Now d(P n , J XY ) = α(P n , J XY ) n -λ 2 (J XY ) = X • Y n - X k • Y k = 0 .
(5.9)

For J Xy we have

X + a Y k , Y + b X k ∈ J Xy =⇒ b X ≤ y .
(5.10)

If we write J y to denote the one-dimensional interval [0, y], we observe that

α(P n , J Xy ) = α (b(X), J y ) (5.11) and λ 2 (J Xy ) = X k • y k = Xλ 1 (J y ) n (5.12)
Therefore: . Note, however, that such a construction does not extend automatically to higher dimensions. We denote P S a point set obtained using this construction.

d(P n , J Xy ) = α(P n , J Xy ) n -λ 2 (J Xy ) = α(b(X), J y ) n - X • λ 1 (J y ) n = X n A(J y ; b(X)) X -λ 1 (J y ) = X n d(b(X), J y ) ≤ X n D * (b(X)) . ( 5

This lemma bounds the discrepancy of our construction and its corollary

Lemma 2 If D * (b (X) ) has a discrepancy in O(log(n)/n) and D * (a (Y ) ) also has a discrepancy in O(log(n)/n), then D * (P n ) is in O(log(n)/n).

Template point set

In practice, we chose the Van der Corput sequence φ (see Section 3.3.2 on page 64). The resulting template LD point set is then defined by:

P φ := (X + φ(Y ), Y + φ(X)) : X, Y ∈ {0 . . . k -1} .
(5.17)

This set, denoted P φ (Fig. 5.3 on the facing page), will serve as our template LD set. Thanks to the stratification, our construction is inherently indexed in two dimensions, using a pair of indices. This allows to align our LD set with any stratified set, on a one-to-one basis. This will be a key element to control the spectral properties of the set.

Following Lemma 2 on the previous page, P φ is a low discrepancy point set with discrepancy O (log(n)/n). It is worth noting that P φ coincides with the Hammersley construction when the generated number of points is a power of two; but not otherwise.

Once we have this initial set, we can optimize it by scrambling the samples to match a given spectral target. The critical part of this optimization is to do this without affecting the discrepancy of the complete set.

Scrambling

We can randomly scramble our template set by applying a discrepancy-preserving rearrangement of the sequences S used as offsets for each row and column. Here, the sequence S is the van der Corput sequence.

Unlike the global bit-wise scrambling techniques in common LD constructions, we randomly permute the entries of the constituent sequences (in rows and columns) locally within small . chunks (Fig. 5.4). We choose a chunk size, m, and re-order every subsequent m entries of the input sequence S (we assume that m divides k): Note that the statement of the mentioned theorem compares corresponding entries in two finite sequences. But since we are referring to subsets, these can be arranged in any order without altering their discrepancy. To give an example, the sets (0, 0.25, 0.5, 0.75) and (0, 0.5, 0.25, 0.75) have different discrepancies as sequences, but the same discrepancy as fixed sets.

S j•m+i =S j•m+σj (i) ; j ∈ {0 . . . k/m} ; i ∈ {0 . . . m -1} , ( 5
To incorporate this 1-D chunk-wise permutation into P φ , we apply two passes: first along the X (columns), then along the Y (rows). For each pass, we divide the X (resp. Y ) into chuks of size m. Then, within each chunk, we reorder the offsets b (resp. a). Note that along the rows we reorder the vertical offsets and along the columns, we reorder the horizontal offsets. Algorithm 3 sums up this whole permutation.

At this point, the resulting point set (Fig.

Blue Noise scrambling

Instead of considering random permutations, we can construct the permutations such that the point set P φ is close to a given reference point set. Consider a stratified t × t point set (t=128 in our experiments), optimized to a specific blue noise profile (BNOT [START_REF] De Goes | Blue noise through optimal transport[END_REF],

Step [Hec+13], .

Poisson Disk, etc). We assume that t is a multiple of m, so that each row and column of the point set can be decomposed into chunks of size m.

The main idea is to start with a t×t chunk of the P φ template set, and to find the permutations for all 1-D chunks that would transform this point set so that it approximates the reference set. For each chunk of size m from the template LD set and the reference set, we rearrange the points of the template set so that the smallest offset a of the template aligns with the smallest offset a of the reference, the second smallest with the second smallest, and so on. This process thus minimizes the distance between the two sets. Figure 5.6 illustrates this for a single 1-D chunk.

The new sequence of indices, is then appropriately encoded and stored in the look up table. For a given point in a stratum (X t , Y t ) in the t × t domain, the look-up table entry LU T (X t , Y t ) contains the offset (X m , Y m ) ∈ m × m, of the new sequence numbers within the m × m block.

The offset of the sample in (X t , Y t ) cannot however simply be replaced by the offset read from LU T (X t , Y t ). Doing this would replicate the t × t permuted point set, which would break the low discrepancy property since several points would project to the same 1-D point on either of the two axes. To overcome this issue, we define an indexing scheme related to the underlying 1-D LD sequence and the bit-reversal principle.

For a given abscissa X, we define X t := X mod t, its position in a t × t domain, and X m := X mod m, its position in a chunk of size m within the t×t domain (we assume similar definitions for Y t and Y m ). The final offsets are given by:

a := φ(Y -Y m + L Y ) (5.19) b := φ(X -X m + L X ) , (5.20)
where (L X , L Y ) := LU T (X t , Y t ).

Eq. (5.19) and (5.20) perform both the initial offset construction, using the van der Corput sequence, and the optimized permutation within chunks of size m. More precisely, for a given X and Y , φ(X) is the vertical offset in the template LD point set (Eq. ( 5.17 on page 98)). The By performing those permutations independently along the x-axis, then along the y-axis, for each m × m block in the point-set, we can create a new pointset P S that is still low discrepancy while presenting the spectral profile of the reference point set. Note that both passes are completely independent, and can therefore be parallelized. Figure 5.8 illustrates this algorithm inside a single m × m block.

5.2 Scrambling X 2 0 2 1 . . . 2 log(n) -1 X mod t { X mod m { X X -X m { L X { φ(X ) 2 -1 2 -2 . . . 2 -log(n)
Note that the parameter m is a trade-off between the degradation of the discrepancy and the desired spectral profile. m = 1 reproduces P S ; m = t produces a merely Latinized copy of the reference set, with limited improvement of its discrepancy. From our experiments, m = 16 is a good compromise. The impact of m on the discrepancy was discussed earlier, following a theorem from [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF] (see Section 4.1.2 on page 75). The size of the reference point set, t 2 , is a trade-off between the memory footprint of the look-up table, and the quality of the spectral profile with respect to the reference set when n > t 2 . Figure 5.11 illustrates how the spectrum evolves with different values for t and m. As can be seen, the bigger the chunks m, the more Blue Noise the spectra, but if m increases too much, the Low Discrepancy property degrades. Also, as we increase t, we need less and less replication to generate our point sets which leads Algorithm 4: Generating an optimized LD set from a reference stratified point set and our template LD point set.

1 Initialize the output set to P S from Eq. (5.17); 2 for each column X do 3 for each chunk j do 4

Retrieve the vector

L x = {a X,Y } j•m+m-1 Y =k•m
from the reference set; 5 Retrieve the corresponding vector, Q x , from the output set;

6

Rearrange Q x so that the smallest entry aligns with the smallest entry in L x , the second smallest with the second smallest, and so on; to cleaner spectra but at the cost of an increase in the memory requirements. In our tests, the parameters m = 16 and t = 128 are a good trade-off.

Using previous notations, we define our final Low-Discrepancy Anti-Aliased sampler as follows:

LDBN := X + φ(Y -Y m + L Y ), Y + φ(X -X m + L X ) . (5.21)
Following Lemma 2 on page 97, we can affirm that the LDBN point set defined from van der Corput sequences in Eq.( 5.21) has a discrepancy in O log(n) n , and is thus a low discrepancy set.

We further empirically demonstrate in Chapter 7 on page 133 that when optimizing the look-up table with a reference point set with specific spectral properties (e.g, Step blue noise), LDBN preserves the low discrepancy of P φ while presenting a different spectral profile. The result of this sampler with different spectral profiles are presented in Figure 5.9 on the preceding page, and in Figure 5.10 on the previous page. Algorithm 4 sums up our optimized permutation.

Conclusion

On the practical side, LDBN is extremely computationally efficient, since it can be implemented using only bit manipulations and memory lookups. It is also considerably efficient in storage, since only integer permutations in [0, m) have to be stored, rather than the actual real-valued offsets (a, b). If m = 16, then the storage space reduces to only one byte per entry.

However, despite being fast and generating Low Discrepancy sets with controllable spectrum, this method is still limited. It only operates on 2-D points, and only generates sets of samples.

Multi-Dimensional Sequences with Low-Discrepancy Blue-Noise 2-D Projections

In this chapter, we detail a new 2-D scrambling that can be applied on projections of the Sobol sequence. This lets us generate s-dimensional sequences of samples with 2-D projections that are both low-discrepancy and present a Fourier spectrum approaching the Blue Noise spectrum.

Our tile-based method performs permutations on subsets of the original Sobol sequence. In the large space of all possible permutations, we choose those which better approach the target Blue Noise property, using pair-correlation statistics. We pre-calculate such optimized permutations for each possible Sobol sub-pattern, and store them in a Look-Up Table (LUT) efficiently accessible at runtime. We provide a complete and rigorous proof that such permutations preserve Sobol's net properties and thus the LDS properties of the point set in 2-D projections. Our muti-dimensional construction is computationally efficient, has relatively low memory footprint and supports adaptive sampling.

We start by presenting a general overview of our permutation (Section 6.1). Then we will dive into deeper details (Section 6.2 on page 117), wich will let us formally demonstrate its low discrepancy and sequentiality (Section 6.3 on page 122).

Overview

In this section, we provide a high-level, simplified and intuitive explanation of our scrambling. Rigorous definitions and technical details are provided in Section 6.2 on page 117. Section 6.1.1 presents our 2-D scrambling, which relies on a precomputed set of optimized local permutations. We then describe how we can build those permutations, along with the performances of the associated Look-Up Table in Section 6.1.2 on page 110 and 6.1.3 on page 114. Finally, we will see how this 2-D scrambling can help creating s-D sampling patterns Section 6.1.4 on page 114.

2-D Scrambling

Our sampler relies on independent scrambling of 2-D projections of the Sobol sequence. In this section we can therefore start by presenting how this scrambling applies on a single pair of dimensions without loss of generality. Our construction can be conceptually considered as a set of staged local permutations of tiles of K 2 samples of the initial Sobol LDS (the leftmost column in Figure 6.1). The factor K = 2 n , n ∈ N is used to determine the number of samples of those tiles. Our basic 2-D scrambling is illustrated in Figure 6.1.

P 0 S 1 T 1 0 { π 0 Q 1 π 1 0 π 1 1 π 1 2 π 1 3 π 1 4 π 1 5 π 1 6 π 1 7 π 1 8 π 1 9 π 1 10 π 1 11 π 1 12 π 1 13 π 1 14 π 1 15 P 1 S 2 π 0 π 1 0 π 1 1 π 1 2 π 1 3 π 1 4 π 1 5 π 1 6 π 1 7 π 1 8 π 1 9 π 1 10 π 1 11 π 1 12 π 1 13 π 1 14 π 1 15 Q 2 P 2 (a) (b) (c) (d) (e) (f) (g)
On the left column of this figure, we present the original Sobol point sets (upper row: 16 1 points; middle row: 16 2 points; lower row: 16 3 points), generated using Sobol polynomials x 2 + x 1 + and x 4 + x 3 + 1 for dimensions x and y. Note that all square partitions delimited by thin blue lines contain exactly 16 points each; this remains true for any power p of 16 p sampling points.

The upper row is the first iteration of our algorithm. We generate P λ for λ = 1. The permutation π 0 makes first 16 points of the Sobol sequence more evenly distributed over the domain [0, 1) 2 .

The second iteration is shown in the middle row. We generate P λ for λ = 2, iterating over P 1 . The permutation π 0 applied to 16 2 points of the Sobol sequence keeps the first 16 points welldistributed (large colored dots in the middle sub-figure, middle row), whereas the remaining points (smaller blue dots) are unevenly distributed. Staged π 1 i permutations make all first points of the the sequence more evenly distributed (rightmost sub-image, middle row). Pivot points (the first point of each group of 16 points, according to Sobol's numbering) are shown as colored larger dots for the first 16 2 points and larger red dots for the first 16 3 points of the sequence. Note that after application of π 0 on 16 2 points, the pivots of permuted tiles of points match the distribution of the first 16 points after application of π 0 .

Finally, the lower row present a thrid iteration. We generate P λ for λ = 3, iterating over P 2 . Similarly, 3 levels of permutations are needed in order to obtain even distribution of first 16 3 points of the sequence. Note that after application of π 0 on 16 3 points, the first 16 points are evenly distributed, as in upper row, whereas staged π 1 i applied on 16 3 points keeps the first 16 2 points evenly distributed, as in the middle row. Only three levels of staged permutations are required to achieve an even distribution of all 16 3 points (rightmost sub-image, lower row). Please zoom in the sub-figures of the lower row to see fine details.

In this illustration, we work with tiles of K 2 = 16 samples. Note that similar constructions exist with tiles of K 2 = 64, K 2 = 256 or any K 2 with K = 2 n , n ∈ N * samples. In this illustration, we also used the pair of indices 3 and 7 for the Sobol sequence but note that that other pairs of indices can be used.

Our construction works by iteratively enriching an initial set. It starts at level λ = 0 with a set of K 2 Sobol samples. For those first K 2 points of the initial set (the top row in Figure 6.1 on the facing page), a good pair of Owen's permutation trees [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF], chosen according to a given blue-noise criterion, can be found (see Section 6.1.2 on the next page for more details on how we compute them). Our pattern is then permuted using Owen's scrambling with those optimized permutation trees, to produce a more even distribution of sampling points (the rightmost subimage of the top row in Figure 6.1 on the facing page). Note that thanks to the properties of Owen's scrambling, this permutation preserves the t factor of the initial Sobol set of K 2 samples. (see Section 6.3 on page 122). In Figure 6.1 on the facing page, we denote this permutation π 0 . Please refer to Section 4.1.3 on page 76 for a detailed description of Owen's scrambling.

We then refine our set of K 2 samples to create a new set of K 2•2 = 256 samples. We do this by starting over from the Sobol set containing K 2•2 = 256 samples (the middle row in Figure 6.1 on the facing page). We first apply the same permutation π 0 on the entire point set. Then, we subdivide our domain into K 2 tiles, containing K 2 samples each (see Section 6.2 on page 117). For each tile, denoted T λ i (λ = 1), we apply a new optimized permutation π 1 i . This generates an even distribution of all K 2•2 samples, as shown in the rightmost sub-images of the middle row in Figure 6.1 on the facing page. Note that after the first stage of permutation π 0 , the K 2 = 16 first samples of the set are the same as the optimized samples from the previous optimized set (shown as large colored dots in the Q 1 and Q 2 images from Figure 6.1 on the preceding page). As the second stage of permutations π 1 i also does not move those first K 2 samples, we will refer to those points as pivots, because they are immovable samples w.r.t. previous stages.

It then goes similarly for the following iteration (the lower row in Figure 6.1 on the facing page). We start with the set of K 2•3 = 4096 Sobol samples. We then undergo 3 stages of permutations, in order to keep first K 2 and K 2•2 samples well-distributed, and achieve good distribution of all K 2•3 samples. This process can be continued to higher orders λ of K 2(λ+1) samples. Note that the first λ stages of permutations (that keeps the distribution of the first K 2 , K 4 , • • • , K 2λ samples) can be applied with a single xoring between the initial Sobol point set from the previous level S λ-1 and the optimized point set P λ-1 obtained from this Sobol set (this will be detailed in Section 6.2 on page 117). Applying those xor vectors on the initial Sobol set of the current level S λ creates a temporary set Q λ ,. This set is then locally scrambled to generate the optimized set P λ . This process will be detailed in section 6.2 on page 117.

Algorithm 5: Refining P λ-1 to P λ input : Point Set P λ-1 at level λ -1 and the lookup table LU T (•). output : Point Set P λ . // Construct Q λ from S λ-1 , S λ , and Algorithm 5 sums up the general framework for the application of our optimized Owen's permutations on Sobol point sets. We will now see how we can carefully choose our permutation trees, to ensure that our 2-D scrambling will have a spectrum close to the Blue Noise spectrum while ensuring that our local scramblings will not change the position of the points that were generated at the previous level (to preserve the sequentiality).

P λ-1 1 forall tiles T λ i do 2 s (i) ← S λ-1 ∩ T λ i ; 3 p (i) ← P λ-1 ∩ T λ i ; 4 v λ r ← s (i) ⊕ p (i) ; 5 forall points s ∈ S λ ∩ T λ i do 6 Q λ ← Q λ ∪ (s ⊕ v λ i );

Local permutations with spectral control

In this section, we see how we can choose our 2-D local permutations to ensure that our scrambling preserves the discrepancy while creating a Blue Noise spectrum. Our permutations are based on Owen's scrambling [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF] (Section 4.1.3 on page 76) which means that a 2-D permutation is represented as a pair of Boolean trees. Our goal is thus to figure how to choose the flags of those trees, in order to both preserve the sequentiality of the whole global set and create locally Blue Noise distributions of samples. To do so, we follow a method similar in essence to the one presented in Section 4.4.2 on page 87 (that optimizes the Boolan trees using genetic algorithms). However, in our particular case of local independant permutations, further constraints are needed to ensure that the discrepancy of our global set is preserved.

Let us first consider that we have a point set Q λ with the following property: each tile T λ i contains K 2 samples (we remind that K = 2 n , n ∈ N), with a single pivot sample (see Figure 6.1 on page 108). We further assume that Q λ is a (t, λ, 2)-net in base K 2 , and that each pattern inside each tile T λ i of Q λ is a (t, 2n, 2)-net in base 2 (we will come back to this in Section 6.2 on page 117 and 6.3 on page 122).

To obtain the point set P λ , we locally scramble the points inside each tile T λ i of Q λ using optimized permutations. Each tile T λ i contains a different set of samples. Each unique set of samples requires its own permutation π λ i To preserve sequentiality and discrepancy of the original 2-D set, each optimized permutation π λ i needs to fulfil the three following conditions (see Lemmas 4 on page 126 and 3 on page 125):

(i) π λ i does not move the pivot sample.

(ii) π λ i is (t, 2n, 2)-net preserving in base 2.

(iii) π λ i does not change the 1-D projections of the input pattern.

Property (i) is mandatory to ensure the sequentiality, Property (ii) is mandatory to ensure that the discrepancy is preserved within each tile, and Property (iii) ensures that the discrepancy of the global 2-D set is also preserved. We say that a permutation π λ i is an admissible permutations for a given pattern of samples if it fulfils the above mentioned constraints.

Our permutation applies an Owen's scrambling over the point set within T λ i . If we consider any set of s static Boolean trees of depth m, Owen's scrambling is known to be (t, m, 2s)-net preserving in base 2. Therefore, any permutation with s = 2 and m = 2n, validates property (ii). However, it does not necessarily preserve the 1-D projections of samples (property (iii)) when coordinates are expressed on more than m digits, or when there are fewer than 2 m samples. To alleviate this, we apply the scrambling on the first m bits and exchange the trailing ones in order to preserve the 1-D projections (see Section 6.2 on page 117). As a consequence, for each tile T λ i , any such tree with m = log(K) 2 induces a permutation π λ r satisfying (ii) and (iii).

To satisfy (i), we need to be able to scramble the point set while preserving the position of the pivot sample. This can be done using Owen's scrambling. If we list all the possible branches within the permutation tree, each of them will apply on a single pair of samples. Therefore, if along the path leading to the pivot sample, we set all the flags to 0, this sample will not be permuted when all others will move. We thus have m flags set to 0, while the remaining K 2 -1m flags are free, which gives us enough degrees of freedom to still have a control over the Fourier spectrum.

Using such admissible permutations, we can generate a multidimensional sequence with low discrepancy 2-D projections presenting a stratified spectrum (Fig. 6.2 on the next page). The computations of such random permutation trees is quite efficient as we simply fill the tree with random values apart from a branch that is fixed to 0. Furthermore, applying such is as computationally efficient as the original Owen's scrambling with finite small permutation trees. However, it does not lead to the Blue Noise spectrum we expect.

We will denote an ad-hoc permutation for a given pattern as an admissible permutation for this pattern that leads to a Fourier spectrum as close as possible to a Blue Noise spectrum.

Identifying such permutations among the set of all admissible permutations of a pattern can be very challenging. We look for transformations that optimize the spectral content of the point set, e.g. by increasing the minimum distance between points or obtaining a pattern who is as close as possible to a blue noise one. To characterize a pattern as a realization of a stationary and ergodic stochastic point process, we consider its pair correlation functions (PCF) which can be related to the spectral content of the point set (Section 2.3.1 on page 26 and [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF]).

To evaluate the quality of a pattern with respect to a blue noise one, we compute the l ∞ or l 2 norm of the two PCFs as discussed in [START_REF] Öztireli | Analysis and Synthesis of Point Distributions based on Pair Correlation[END_REF] (in our experiments, both measures give similar results).

In dimension 2 for a given K, the point set {q} λ i (defined as the sub set of samples from Q λ that belong to a tile T λ i , formally

{q} λ i = Q λ ∩ T λ i ) contains K 2 sample.
For each tile, the permutation is encoded with s Boolean trees of depth m with (K 2 -1m) free flags. This gives us a total of 2 K 2 -1-m permutations among which we want to find the ones that lead to the best Blue Noise spectrum.

As an example, for K = 4, we have 2048 2 different admissible permutations to explore for a single pattern. This lets us perform an exhaustive listing of all those permutations on the GPU and keep the best one (which takes about 1 second on a GeForce GTX 760). For each permutation, we evaluate its quality using its PCF. We compare its PCF to the PCF of Blue Noise, and measure the l ∞ distance between the two. The one with the lowest distance is kept.

For K = 4 (16 points), an exhaustive search is possible. For K > 4, we optimize exhaustively the first 16 bits of the permutation tree and fill the underlying digits at random. As our GPU implementation performs this search in 1s, we can brute force the last digits to ensure that, even when randomly chosen, they do not harm the quality of the pattern. This whole process takes up to 2s, still on a GPU GeForce GTX 760. Figure 6.3 and 6.4 present the characteristics of our scrambling when using ad-hoc permutations (for K = 4 and K = 8). Note that applying a given permutation is still very efficient. However, computing this permutation has become quite costly. Especially if you point out that generating K 2n samples implies that we have 1 +

K 2 + K 4 + • • • + K 2(n-1
) tiles to scramble. However, as will be detailed in Section 6.1.3 on the next page, the recursive structure of Sobol implies that several tiles present similar sets of samples, allowing to reuse at times previously optimized permutations. A smaller K value increases the amount of similarity within tiles.

Once we have computed the set of ad-hoc permutations regarding a specific pattern, we store them in the lookup table, indexed by the pattern {q} λ i . Note that multiple ad-hoc permutations can be stored for the pattern {q} λ i . This allows us to randomize the scrambler by randomly picking a good permutation (see Algorithm 5 on page 110). The number of entries in the LUT is the number of distinct {q} λ i patterns.

Look Up Table

The LUT stores the ad-hoc permutations for each pattern {q} λ i in order to generate the set P λ . Since the set Q λ is defined from P λ-1 (Fig. 6.1 on page 108), we fill the LUT while following the hierarchical construction. We start from λ = 1, and we generate successive P λ and Q λ sets. Each time we have to permute a pattern {q} λ i , whose optimized permutation is not in the LUT, we explore admissible permutations and add the best one to the LUT (see Section 6.1.2 to see how we pick this best one).

All patterns {q} λ

i have the same number of points K 2 . The total number of possible distinct patterns, (and thus the upper bound on the LUT size) is difficult to estimate (e.g. O((K 2 !) 2 ) if we enumerate the set of Latin hypercube configurations, from which we must remove non-dyadic sets). However, as we start from the Sobol sequence which present a very repetitive structure, the number of distinct sets is limited. Empirically, as λ increases, the number of entries in the lookup table stabilizes quickly. As shown in Figure 6.5 on the facing page, for K = 2 only 12 configurations exist up to 10 7 samples (λ = 10). For K = 4, only 512 configurations exist up to 10 6 samples (λ = 6). For K = 8, we found up to 1899 configurations for about 2.10 5 samples (λ = 3), but the table does not seem to have reached saturation (Fig. 6.5 on the next page)).

In 2-D, a permutation is defined by 2 Boolean trees with K 2 -1 nodes for K := 2 n . Hence, for K = 4 (resp. K = 8), a straightforward encoding of a permutation π λ i requires 30 bits (resp. 126 bits). Each entry in the lookup table is indexed by a pattern containing K 2 samples, each sample being encoded onto 2 • 2 • n bits. If we allow the LUT to contain m ad-hoc permutations for each pattern, we need

K 2 • 4 • n + m • 2(K 2 -1)
bits.

For K = 4 and m = 1, the LUT size is 512 • (86 + 30) = 59392 bits which is about 7 kbytes (about 485 kbytes for K = 8).

Note that this LUT contains optimized permutations for 2-D sampling patterns only. Unfortunately, our scrambling does not properly extend to higher dimensions (more details in Section 6.2). Therefore, to use our method with s-D point set, we can only use our optimized scrambling in some 2-D projections.

Using Our Sampler in Higher Dimensions

In many rendering applications, higher dimensional samples may be required. The permutation described above is only defined in 2-D. However, even though it only performs operations that are perfectly defined in s-D, our construction cannot be straightforwardly extended to s-D, as it would no longer preserve the discrepancy. The condition (iii) for admissible permutations only holds for 2-D samples. We instead extend our method to s-dimensions by relying on independent scrambling of several 2-D projections (see Figure 6.6 on page 116). In short, we can use a s-D Sobol sequence, optimize some 2-D projections with respect to pre-computed lookup tables and perform an on-the-fly random permutation of the remaining dimensions. Such an approach is made possible since the initial sequence is consistent through the dimensions and has LD properties. The selected 2-D projections will have low discrepancy properties and an optimized spectra. In rendering applications, selected projections could be the image plane, the lens, area lights, etc.

In order to have the (t, λ, 2)-net in base K 2 property on the optimized pair of dimensions, we need each set {s} λ r := S λ ∪ T λ i to be a (t, ns, 2)-net in base 2 in these dimensions (see Prop. 4 on page 126). In our implementation, we have used certain combinations of primitive polynomials' indices, found empirically. When K=4, we used for our 4 dimensions the indices (1, 2) and the indices (3,7). When K=8, we used for our 6 dimensions the indices the indices (1, 4), [START_REF]D point set with stratified 1-D projections[END_REF]3) and the indices (5,7). Note that as K increases, more and more dimensions can be optimized.

For the remaining dimensions, when the number of samples is fixed, one can use a classical Owen's scrambling of the point coordinates with a tree of fixed depth (Fig. 6.7 on the following page). In our hierarchical setting when the number of points is unknown and increases on the fly, we simulate a tree of infinite depth by iteratively expanding a tree of fixed depth, as detailed in Section 4.2 on page 79.

From here, you should have all the intuitions required to understand how our method works and why it works. However, to formally prove how our scrambling preserves the Low Discrepancy property of Sobol projections, we need a much more formal and detailed description of our sampler. This will be the focus of the next section. We will further validate our method in Chapter 7 on page 133 by performing spectral/discrepancy/aliasing analysis of the achieved distributions, and provide variance analysis for several target integrands of theoretical and practical interest.

Technical description

In this section, we will present our method in much deeper details. Thanks to those formal definitions, we will be able to demonstrate that it preserves the discrepancy of the samples in Lemmas 3 on page 125 , 4 on page 126 and 5 on page 129 (Section 6.3 on page 122). Note that we will detail our method in a supposedly s-D context, but in practice, it only preserves the discrepancy for 2-D point sets. We also stress that all illustrations in Section 6.2 and 6.3 are given for K = 2 (4 samples per tile), to ease understanding. In practice, our method is usually applied with K = 4 and K = 8.

As a short preliminary, we remind that we denote the binary representation of a sample as follows. Any s-D sample x ∈ [0, 1) s can be expressed as (x 1 , x 2 , . . . , x s ) with

x d := ∞ i=1 x i d 2 -1-i x d = .x 1 d x 2 d x 3 d . . . ,
where x i d is the i th bit of x d . Note that the samples are defined in [0, 1) s , therefore, x 1 d is the most significant bit, meaning x 1 d corresponds to 0.5, x 2 d to 0.25, . . .

We stress that the proofs in this section require knowledge of the notations and definitions of (t, k, s)-nets and sequences that were presented in Section 3.3 on page 57.

Domain subdivision

Our scrambling considers a regular subdivision of the domain [0, 1) s with a subdivision factor K := 2 n , n ∈ N * . We denote λ the level of subdivision of the domain (Fig 6 .8 on the next page).

This subdivision leads to square tiles. Each time we subdivide the domain, we subdivide the tiles at level λ by the factor K, creating K (λ+1)s new tiles at level λ + 1. A tile T λ r is a s-D interval of space defined at a level λ spanning over

r d K λ , r d +1 K λ s
. The union of all the tiles at a level λ forms a partition of the domain [0, 1) D (Fig 6.9 on the following page).

To sample the domain subdivided at a level λ, we generate the first K s(λ+1) = 2 ns(λ+1) samples from the Sobol sequence [START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] (Fig 6 .10 on the next page). By definition of the Sobol sequence, S λ is a (t, k, s)-net in base b = 2 (Section 3.3 on page 57). Since S λ contains 2 ns(λ+1) samples, we have b k = 2 ns(λ+1) and therefore S λ is a (t, ns(λ + 1), s)net in base b = 2. Our set S λ (defined as the set containing the first K s(λ+1) samples of the sobol sequence with K := 2 n ) is a (t, ns(λ + 1), s)-net in base b = 2. In our case, we assume that the value t of our (t, k, s)-nets is the minimal possible value, for which (and over which) all intervals are guaranteed, we have the following property:

K = 2, λ = 0 K = 2, λ = 1 K = 2, λ = 2
Property 1 ∀λ ∈ N, if S λ is a (t, ns(λ + 1), s)-net in base b = 2, S λ is a ( t ns , λ + 1, s)-net in base K s , with K := 2 n .
From Property 1, we have that if the t value for the Sobol sequence verifies t ≤ ns, we are guaranteed to have K s samples inside each tile T λ r . If t > ns, we can still have K s samples in each tile, there is just no guarantee for it. Note that in practice, we also never have t ≤ ns.

S λ-1 T λ-1 r { r {s} λ-1,λ r T λ r { S λ T λ r { r {s} λ r T λ r {
Fig. 6.11.: s = 2, K = 2, λ = 2. {s} λ r is the intersection between S λ and the tile T λ r . Similarly, {s} λ-1,λ r is the intersection between S λ-1 and the tile T λ r . We now define {s} λ r := S λ ∩ T λ r the set of samples from S λ that belong to the tile T λ r . We also define {s} λ0,λ1 r := S λ0 ∩ T λ1 r the set of sample from S λ0 that belong to the tile T λ1 r (see Figure 6.11).

It is important to note that since the samples of S λ come from the Sobol sequence, they are indexed [START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. We will denote s (i) the i th sample of the sequence S λ . A reader familiar with the Sobol sequence may note that the samples of {s} λ r are also indexed but their indices i are not consecutive (Fig 6 .12).

We then permute the samples from S λ , where S λ is a set containing the first K s(λ+1) samples of the Sobol sequence. We do this in such a way that the newly generated set of samples P λ in [0, 1) s has the following properties:

• If s = 2, ∀λ, if S λ is a (t, λ + 1, s)-net in base K s , P λ is a (t, λ + 1, s)-net in base K s
(Lemma 5 on page 129).

• ∀λ, P λ presents a user defined Fourier spectrum.

• ∀λ > 0, P λ-1 ⊂ P λ (otherwise P wouldn't be a sequence), and the samples in P λ are indexed (Lemma 6 on page 129).

To ensure those properties, our algorithm creates a set P λ by permuting the set S λ while taking as input the two sets, P λ-1 and S λ-1 . Note that we define P -1 as the point set containing a single sample s = 0, and S -1 = P -1 . The sets S λ and S λ-1 are coming from the same specific Sobol sequence, with indices chosen such that S λ will be stratified ∀λ. The use of the Sobol sequence also ensures that S λ-1 ⊂ S λ .

Then, our permutation is done in two steps (Fig. 6.1 on page 108)

1. A first permutation, denoted Ψ, applies local digital shifts on S λ to create a new point set

Q λ such that P λ-1 ⊂ Q λ .
2. A second permutation, denoted Π, applies local permutations π λ r on Q λ to create a new point set P λ such that it presents the targeted Fourier spectrum.

Creating

Q λ from S λ and P λ-1
The aim of our first permutation Ψ is to create a set Q λ from S λ , so that P λ-1 ⊂ Q λ . To achieve this, we first compute a set of vectors, such that each vector represent the displacement between each sample of S λ-1 and P λ-1 . Since the samples are indexed, both in S λ and in P λ , we can match those sets on a one to one basis. This will let us xor each sample of S λ-1 with its corresponding sample in P λ-1 (the one with the same index). Therefore, our new vector set, denoted V λ-1 is also ordered, and its i th vector, denoted v (i) is formally defined as v (i) := s (i) ⊕ p (i) , (

where s (i) is the i th sample of S λ-1 and p (i) is the i th sample of P λ-1 .

S λ-1 is stratified with a single sample inside each tile T λ r . Thus, each vector v (i) is associated with one and only one tile T λ r . As each vector v (i) is composed from a sample s (i) , which uniquely belongs to a tile T λ r . We denote v λ r the vector v (i) associated with the tile T λ r . We can now define formally V λ-1 as

V λ-1 := r v λ r , (6.2)
which is illustrated Figure 6.13 on the next page.

This leads to the formal definition of Q λ as where we define the ⊕ operation between a set of sample X and a single vector v as X ⊕ v := {x (i) ⊕ v} where x ∈ X . This is illustrated Figure 6.14. We also note that each pattern {s} λ r is xored with a single vector to create {q} λ r , thus, as the xor operator is net preserving [START_REF] Lemieux | Monte Carlo and Quasi Monte Carlo Sampling[END_REF], we have the following property for

Q λ := r {{s} λ r ⊕ v λ r }, (6.3)
S λ , P λ-1 V λ-1 Q λ , P λ-1
Q λ Property 2 If ∀T λ r , N ({s} λ r ), is a (t, ns, s)-net in base 2, then ∀T λ r , N ({q} λ r ), is a (t, ns, s)-net in base 2.
Then, from the point set Q λ , we can create our optimized point set P λ by locally scrambling the samples.

Creating

P λ from Q λ
Then, we permute the samples inside each tile T λ r of Q λ to create P λ . Note that there are K s samples of Q λ inside T λ r , forming a point set denoted {q} λ r . Each set {q} λ r contains a unique sample p such that p ∈ P λ-1 . Note that each set {q} λ r contains samples with coordinates that only cover the tile T λ r (Fig. 6.15 on the next page). We denote N ({q} λ r ) the operation that scale this set to the domain [0, 1) s , with

N -1 (N ({q} λ r )) = {q} λ r . {q} λ r N ({q} λ r ) π λ r (N ({q} λ r )) {p} λ r := N -1 (π λ r (N ({q} λ r )))
Fig. 6.15.: s = 2, K = 2, λ = 2. We permute each tile T λ r from Q λ to create a new arrangement of samples within each tile with the desired spectral properties.

Once we have this set, scaled to the domain [0, 1) s , we can associate with each tile T λ r a permutation π λ r , that have the following properties (denoted as (i), (ii) and (iii) in Section 6.1):

Property 3 If N ({q} λ r ), is a (t, k, s)-net in base 2, N (π λ r ({q} λ r )) is also a (t, k, s)-net in base 2. Property 4 If p ∈ {q} λ r and p ∈ P λ-1 , p ∈ π λ r ({q} λ r ).
Property 5 π λ r preserves the 1-D projections of {q} λ r .

We can now formally define P λ with

P λ := r N -1 (π λ r (N ({q} λ r ))) (6.4)
A possible construction algorithm to generate admissible permutations π is given in Section 6.1 on page 107 and is detailed in Section 6.3.6 on page 130.

Proofs

We will now prove that our scrambling preserves the net properties of S λ . We will first define some notations, and then we will prove that both Ψ and Π are (t, λ+1, 2)-net preserving in base K s . Finally, we will prove that P λ does contain the samples of P λ-1 , identically indexed.

Notations

In this section we present all the notations that we use in the following proofs. Please note that even though our permutation is valid in s-D, it only preserves the net properties in 2-D. Therefore all those notations consider a 2-D domain.

Each 2-D sample (x, y) in a point set X , where all the samples in X are defined on ns(λ + 1) bits, and where X is a (t, λ + 1, 2)-net, can be expressed in binary as

x = ns(λ+1) i=0 x i 2 -1-i , y = ns(λ+1) i=0 y i 2 -1-i ,
where x i is the i th bit of x and y i is the i th bit of y.

x = ns(λ + 1)

x 0 { x i { y = ns(λ + 1) y 0 { y i { Fig. 6.16.: Binary representation of a sample (x, y) ∈ [0, 1) s
We define x| a (resp. y| a ) as a binary truncating of the x (resp. y) coordinate of a sample (x, y) where

x| a = a i=0 x i 2 -1-i , x| a = ns(λ + 1)
x a { Fig. 6.17.: Binary representation of x| a

We extend this definition to (x, y)| a , for (x, y) ∈ X in a specific way. We thus define (x, y)| a as (x, y)| a = (x| a , y| ns(λ+1-t)-a )

This means that we truncate a total of ns(λ + 1t) bits from (x, y), truncating a bits from x and ns(λ + 1t)a bits from y.

Finally, we define X | a as {(x, y)| a , ∀(x, y) ∈ X }, meaning that X | a contains truncated samples from X . Each elementary interval J q0 of X , is of size

1 K sq 0 , 1 K s(λ+1-t-q 0 ) , with 0 ≤ q 0 ≤ λ + 1 -t, q 0 + q 1 = λ + 1 -t ⇔ q 1 = λ + 1 -t -q 0 .
It can thus be encoded in binary with nsq 0 bits for x and ns(λ+1-t-q 0 ) = ns(λ+1-t)-nsq 0 for y. Therefore, for all sample (x, y) ∈ X , their truncated version (x, y)| nsq0 is equivalent to the minimal corner of the elementary interval J nsq0 that contains it.

Following this, a point set

X is a (t, k, 2)-net in base K s if each sample (x, y)| nsi is present K st times in the set X | nsi , with 0 ≤ i < k -t.
(

, y)|nsq 0 q 0 = 0 q 0 = 1 Fig. 6.19.:

K = 2, n = 1. Q 0 is a (0, 1, s)-net in base 2 2 iff ∀q 0 ∈ N with 0 ≤ q 0 ≤ 1, each sample (x, y)| 2 * 2 * q0 , (x, y) ∈ Q 0 appears exactly one time in Q λ | 2•2•q0 .
The direct corollary is that any permutation

Π is (t, k, 2)-net preserving if ∀a ∈ N, a ≤ k -t, Π(X )| a = X | a . ( 6.5)
Finally, we have one last major property, illustrated Figure 6.20 on the facing page.

Property 6 ∀(x, y) ∈ [0, 1) 2 , if (x| a , y| a ) = (x | a , y | a ) then ∀b ≤ a and ∀b ≤ a , (x| b , y| b ) = (x | b , y | b )
This property means that if we truncate the x coordinate (resp. y) of the samples of a set to a (resp. b) digits, and if this truncated set is equal to the truncated set from another point set, then the subsets of those sets that are more truncated, with a < a (resp. b < b), are also equal.

6.3.2 Our permutation Π is (t, λ + 1, 2) preserving in base K 2 Lemma 3 For s = 2, if Q λ is a (t, λ + 1, 2)-net in base K 2 , P λ is a (t, λ + 1, 2)-net in base K 2 .

Proof 2

The point set P λ is defined as Π(Q λ ), where Π applies a scrambling π λ r over each set N ({q} λ r ) from Q λ (Equation 6.4 on page 122). The operator N removes the first nλ bits of the samples in {q} λ r . Therefore, those first nλ bits will be unaffected by the permutation. We apply N -1 by re-adding the first nλ bits of {q} λ r .

In our particular case, the permutations π λ r are not any random permutation but are part of a set of admissible permutations (see Section 6.1.2 on page 110 for details). Each permutation π λ r ∈ Π has two important properties.

First, π λ r preserves the net properties of N ({q} λ r ), and second, each set N ({q} λ r ) is a (t , ns, 2)net in base 2 as it comes from the Sobol sequence. From here, we can translate property 3 on page 122 in binary as

(x| nλ+q 0 , y| nλ+ns-t -q 0 ) = Π (x| nλ+q 0 , y| nλ+ns-t -q 0 ) (6.6) ∀q 0 , q 0 ∈ N, 0 ≤ q 0 ≤ ns -t , q 0 + q 1 = ns -t ⇔ q 1 = ns -t -q 0 .
Furthermore, since π λ r preserves the 1-D projections of samples we have for q 0 = 0, (x| ns(λ+1) , y| nλ+ns-t -q 0 ) = π λ r (x| ns(λ+1) , y| nλ+ns-t -q 0 ) (6. 7)

and (x| nλ , y| ns(λ+1) ) = π λ r (x| nλ , y| ns(λ+1) ) . (6.8)
In 2-D, the elementary intervals in Q λ for q 0 + q 1 = λ + 1t are of size 1 2 nsq 0 × 1 2 nsq 1 with q 1 = λ + 1tq 0 . This means that each elementary interval can be encoded using nsq 0 bits for x and ns(λ + 1tq 0 ) bits for y.

As Q λ is a (t, λ + 1, s)-net in base K s , each sample (x, y)| nsq0 , (x, y) ∈ Q λ appears exactly K st times in Q λ | nsq0 . Therefore, ∀q 0 ∈ N with 0 ≤ q 0 ≤ λ + 1 -t, each sample from the set Q λ | nsq0
stands for an elementary interval.

Proofs

From here, we have 3 possibilities • Either nsq 0 ≤ nλ, in which case from 6 on page 124 and equation 6.8 on the previous page, we have

Q λ | nsq0 = P λ | nsq0
• Either ns(λ + 1tq 0 ) ≤ nλ, in which case from 6 on page 124 and equation 6.7 on the previous page, we also have

Q λ | nsq0 = P λ | nsq0
• Or nsq 0 ≥ nλ and ns(λ + 1tq 0 ) ≥ nλ. In which case, we seek if there is a q 0 such that

Q λ | nsq0 = Q λ | nλ+q 0 , meaning such that nsq 0 = nλ + q 0 and ns(λ + 1 -t -q 0 ) = n(λ + s) -t -q 0 nsq 0 = nλ + q 0 ⇔q 0 = n(sq 0 -λ) then n(λ + s) -t -q 0 = n(λ + s) -t -n(sq 0 -λ) = n(λ + s -sq 0 + λ) -t = n(2λ + s -sq 0 ) -t = ns(λ + 1 - t ns -q 0 )
which means that for t = nst, from equation 6.6 on the preceding page and property 6 on page 124,

Q λ | nsq0 = P λ | nsq0 . Note that as t is the minimal value such that N ({q} λ r ) is a (t , ns, 2)-net, if t ≤ nst, we have that N ({q} λ r ) is a (t , ns, 2) with t = nst.
Therefore from Eq ( 6.5 on page 124), Π preserves the net properties of

Q λ and thus if Q λ is a (t, λ + 1, 2)-net in base K s , Π(Q λ ) = P λ is a (t, λ + 1, 2)-net in base K 2 . 6.3.3 The permutation Ψ is (t, λ + 1, 2) preserving in base K 2 Lemma 4 For s = 2, if S λ is a (t, λ + 1, 2)-net in base K 2 , Q λ is a (t, λ + 1, 2)-net in base K 2 . Proof 3 Q λ is a (t, λ + 1, s) if ∀(x, y) ∈ S λ and ∀(x , y ) ∈ Q λ we have (x| nsa0 , y| ns(λ+1-t-q0) ) = (x | nsa0 , y | ns(λ+1-t-q0) ), (6.9
)

∀q 0 ∈ N and 0 ≤ q 0 ≤ λ + 1.
From Eq ( 6.3 on page 120), we have

Q λ := r ({s} λ r ⊕ v λ r ) (6.10)
For simplicity, we will extend the xor operator ⊕ to indexed sets, by defining that for two indexed sets P, P , defined at a subdivision level λ, with p ∈ P , and p ∈ P with Card(P ) = Card(P ),

P ⊕ P = p (i) ⊕ p (i)
where p (i) and p (i) belong to the same tile T λ r .

Furthermore, we will define the operator φ λ λ . This operator takes as input a stratified point set at a level λ ≤ λ, and duplicates the vector inside each tile T λ r , in order to repeat this vector inside all the subtiles T λ r of T λ r .

More formally, we define φ λ λ so that for two ordered sets, S λ , S λ , where Card(S λ ) = K xs Card(S λ ), x ∈ N, and with s ∈ S λ , s ∈ S λ-1 ,

S λ ⊕ φ λ (S λ-1 ) = r ({s} λ r ⊕ {s } λ r ) (6.11)
where as S λ-1 is stratified at level λ -1, {s } λ r is a singleton.

This operator will let us denote easily the xoring between two point sets defined at different subdivision level (thus with different cardinalities).

From here, we can rewrite Equation 6.10 on the facing page as

Q λ := S λ ⊕ φ λ λ-1 V λ-1 , (6.12)
since the set S λ contains K s(λ+1) samples. From Equation 6.2 on page 120, the set V λ-1 contains K sλ elements, one for each tile T λ r .

With those new definitions, we denote

π λ-1 = V λ-1 ⊕ φ(V λ-2.125 ),
we can thus derive from V λ-1 :

π λ-1 = V λ-1 ⊕ φ λ-1 λ-2.125 (V λ-2.125 ) ⇔V λ-1 = π λ-1 ⊕ φ λ-1 λ-2.125 (V λ-2.125 ) ⇔V λ-1 = π λ-1 ⊕ φ λ-1 λ-2.125 (π λ-2.125 ) ⊕ φ λ-1 λ-3 (V λ-3 ) ⇔V λ-1 = ⊕ λ-1 i=0 φ λ-1 i (π i )
with V -1 = 0.

Proofs

Therefore, Q λ is made from the set S λ through successive ⊕ with φ(π i ). From here, we derive ∀λ , 0 < λ ≤ λ,

S λ -1 ⊕ P λ -1 = V λ -1 ⇔S λ -1 ⊕ φ λ -1 λ -2.125 (V λ -2.125 ) ⊕ P λ -1 = V λ -1 ⊕ φ λ -1 λ -2.125 (V λ -2.125 ) ⇔Q λ -1 ⊕ P λ -1 = π λ -1 ⇔P λ -1 = π λ -1 ⊕ Q λ -1 ⇔Π(Q λ -1 ) = π λ -1 ⊕ Q λ -1
Therefore, xoring Q λ with π λ is similar to applying the permutation Π that created

P λ from Q λ . From Lemma 3 on page 125, if (x, y) ∈ S λ and (x , y ) ∈ S λ ⊕ φ(π λ ), (x| nsq 0 , y| ns(λ +1-t -q 0 ) ) = (x | nsq 0 , y | ns(λ +1-t -q 0 ) ) (6.13)
Furthermore, since π λ = Q λ ⊕ φ(P λ ), when written in binary we have for p ∈ π λ -1

p d = ns(λ +1) i=0 p i d 2 -1-i
where p d is the d th coordinate of p and p i d is the i th bit of p d . We note that ∀i > ns(λ ), p i d = 0.

Geometrically, xoring with π λ means that we move all samples that belong to an elementary interval J q 0 to another J q 0 . However, since all digits of p ∈ π λ are 0 when i > ns(λ ), all samples in intervals that are smaller than 1 2 ns(λ ) on every dimension, meaning all samples from intervals such that nsq 0 > nsq 0 and ns(λ + 1q 0 ) > ns(λ + 1q 0 ), are preserved.

(x| nsq0 , y| ns(λ+1-t-q0) ) = (x | nsq0 , y | ns(λ+1-t-q0) ) (6.14)
Furthermore, since we xor with 0 for every digit > ns(λ ), similarly to the previous proof, we can always reduce the size of the interval on a single dimension, leading to the following equations

(x| nsq 0 , y| ns(λ+1) ) = (x | nsq 0 , y | ns(λ+1) ) (6.15) (x| ns(λ+1) , y| ns(λ +1-t -q 0 ) ) = (x | ns(λ+1) , y | ns(λ +1-t -q 0 ) ) (6.16) ∀q 0 ∈ N, 0 ≤ q 0 ≤ ns(lvl + 1) -n(λ + 1).
From here we have several possibilities for q 0

• q 0 = q 0 , then equation 6.13 equals equation 6.9 on page 126.

• nsq 0 > nsq 0 and ns(λ + 1tq 0 ) > ns(λ + 1tq 0 ) then from property 6 on page 124, equation 6.9 on page 126 is implied by 6.13 on the preceding page.

• nsq 0 < nsq 0 and ns(λ + 1tq 0 ) < ns(λ + 1tq 0 ) then 6.14 on the facing page applies.

• nsq 0 > nsq 0 and ns(λ + 1tq 0 ) < ns(λ + 1tq 0 ) (or the opposite), equation 6.15 on the preceding page (or 6.16 on the facing page) can be used jointly with Property 6 on page 124 to imply 6.9 on page 126

6.3.4 P λ is a (t, λ + 1, s)-net in base K s Lemma 5 For s = 2, if S λ is a (t, λ + 1, 2)-net in base K 2 , P λ is a (t, λ + 1, 2)-net in base K 2 .
Proof 4 P λ is made from S λ by successively applying two permutations. First, the Ψ permutation is applied, creating a point set Q λ , which is a (t, λ + 1, 2)-net in base K 2 from Lemma 4 on page 126. Then, we create P λ by applying Π to Q λ , and as from Lemma 3 on page 125, Π preserves the net properties of Q λ , P λ is a (t, λ + 1, 2)-net in base K 2 .

6.3.5

The samples from P λ-1 are also into P λ Lemma 6 P λ-1 ⊂ P λ Proof 5 Since S λ-1 is a (t, λ + 1, s)-net in base K D from property 1 on page 118 with t < ns and both Q λ-1 and P λ-1 are (t, λ + 1, s)-nets in base K D from Lemmas 4 and 3, for t = 0, each of those sets contains one and only one sample inside each tile T λ r .

S λ-1 = r {s} λ-1,λ r (6.17) Q λ-1 = r {q} λ-1,λ r (6.18) P λ-1 = r {p} λ-1,λ r (6.19)
For the same reason, we have that S λ contains K s samples inside each tile T λ r . As S λ-1 ⊂ S λ from [START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], we have that {s} λ-1,λ r ∈ S λ .

Proofs

A B C D

Initial set Our permutation

When we remove the trailing bits of a sample, i.e. when we quantize this sample to the corner of its elementary interval T λ r (gray cells in above figure), we do not change the LD net properties at level λ. We then apply Owen's scrambling on this quantized set, which is also net-preserving [START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF], Proposition 1). And finally, when we reapply the trailing bits, even if we exchange those bits between the samples, those bits are not significant enough to take a sample out of its elementary interval. As an example, if Owen's scrambling at level λ swaps rows containing A :=(x,y) and B :=(x,y), we construct the points C :=(x,y) and D :=(x,y). This exchange of trailing bits is then net-preserving, and therefore our permutation is net-preserving.

Finally, if the initial set is a (t, k, s)-net, when we quantize it on k digits, all combinations of bits are present. Since Owen's scrambling is net-preserving, after scrambling, all combination of bits will still be represented. Then, simply by reapplying the trailing bits with the proper initial bits, we are guaranteed to preserve the 1-D projections of samples. So our set {π λ r } is a subset of the admissible set of permutation.

Conclusion

Following those proofs, we can fairly say that our new method scrambles Sobol sets to generate 2-D projections with proven Blue Noise and Low discrepancy property. By independently scrambling projections, we also allow the generation of s-D point sets with improved discrepancy and spectrum, as shown in Chapter 7 on page 133. Furthermore, our method is sequential, which would let us iteratively refine the approximation of an integrand if needed.

This projection-wise scrambling is not so much of a drawback for rendering as we solve a recursive set of independent integrands. However, only some particular pairs of Sobol indices can be scrambled with our method, which limits its usage to low dimensional sets (s < 10).

Note that this can be alleviated with our sequential Owen' method (Section 4.2 on page 79), but the spectrum will present a lower quality. Another limitation of our method (which is inherent to the combination of Blue Noise and Low Discrepancy), is that it presents a proper spectrum and discrepancy only for sets of size K 2n , n ∈ N (this is further discussed in Chapter 7 on page 133).

Still, despite the amount of possible future work, our scrambling is the first known solution to create sequential sampling patterns with low discrepancy and improved spectrum in higher dimensions.

Experiments

In this section, we present the results of all experiments we performed when characterizing samplers and validating our methods. All of our evaluations were performed on uniform sampling patterns, the settings for each particular experiments are detailed within the sub sections.

We start by presenting graphs measuring the behaviour of the discrepancy and of the integration error. Those measures are presented for our methods and are thoroughly compared with the state of the art. Then, we present reconstructions and renderings obtained from our methods and some state of the art methods. Unfortunately, due to the complexity of modern rendering systems, this analysis is less thorough but is still very informative. Our final section presents a table of timings for the main samplers.

All those experiments were computed using tools from the UTK Framework [Utk]. For all state of the art samplers, we used, when available, their implementation within the UTK framework.

For other samplers, we used code provided by the authors. Note that for all non-trivial samplers, their implementation within the UTK Framework is a copy of the authors code that was merely adapted to the UTK API. We compare those previously existing samplers to the new methods proposed in Chapter 5 on page 93 and Chapter 6 on page 107, also implemented in the UTK API.

Discrepancy

In this section, we compare the discrepancy of many samplers. Our main interest is to measure how their discrepancy evolves as the number of sample increases, the faster it converges to 0, the better. We start by presenting the l ∞ star discrepancy measure, for 2-D samplers. Then, we will present the l 2 star discrepancy for the same 2-D samplers, and for some 4-D sampling patterns.

l ∞ star discrepancy

In this subsection we present the l ∞ star discrepancy computed using the exact O(n s ) algorithm described in Section 2.1.1 on page 17. The following graphs present the evolution of this discrepancy for many samplers. We first see how Low discrepancy samplers behave (Fig. 7.1 on page 135), then similarly for Blue Noise samplers (Fig. 7.2 on page 135), and finally for samplers that are neither Blue Noise nor Low Discrepancy (Fig. 7.3 on page 136). Then, the discrepancies of the main samplers will be compared to the discrepancy of our new methods (Fig. 7.4 on page 136). Note that all those measures were done with 2-D samplers. On the following graphs, we also plotted the theoretical discrepancies for the 2-D Whitenoise sampler, O(1/n), and for 2-D Low Discrepancy samplers O(log(n)/n).

The first graph (Fig. 7.1 on the next page) presents the discrepancy of various Low Discrepancy samplers. As expected, all curves show the same convergence speed, following the theoretical curve in O(log(n)/n). It is interesting to note how the discrepancy of the Faure [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF] and Sobol [START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] samplers increases and decreases periodically. This is due to the fact that those sequences generate digital nets in base 2. Thus, any number of samples that is not a power of two cannot fill properly the net as there are not enough samples with respect to the number of congruent boxes (see Section 3.3.2 on page 59). For all samplers defined in a base b, we denote octaves the sets presenting b k samples, k ∈ N. Note that in base 2, the Faure and Niederreiter sequences [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] are similar to the Sobol sequence (with a different ordering of samples). Therefore, we plotted the Neiderreiter sequence in base b = 3 here. Also note that in further graphs, we may no longer show the discrepancy outside the octaves. Now if you look at the l ∞ star discrepancy of Blue Noise samplers (Fig. 7.2 on the next page), it can be seen that for many such samplers, the discrepancy is close to the Whitenoise theoretical curve. Note that for such samplers, as they are stochastic, the discrepancy was presented a the maximal discrepancy over around 20 realisations. However, for some methods such as BNOT [START_REF] De Goes | Blue noise through optimal transport[END_REF], CVT [START_REF] Du | Centroidal Voronoi tessellations: Applications and algorithms[END_REF] or CapCVT [START_REF] Chen | Variational blue noise sampling[END_REF], the computational time is really expensive and, for higher number of samples, we sometimes are not able to generate more than one or two sets, hence the increase in variance for the discrepancy measures. In the case of CVT and CapCVT, it can be seen that passed 10000 samples, the discrepancy seems to stop decreasing as the number of samples increases. At this time, we are still not sure to why we observe this behaviour. There are two possible explanations to this phenomenon, either the implementation we have for such samplers (provided by the authors) is slightly flawed, or we are not configuring it properly. If this is not the case, this could be also explained by a unperceived shifting of the points, leading to visually perfect sets, but that is sufficient to affect the discrepancy (similarly to what happens with the sampler from [START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF], as mentioned in Section 3.2.2 on page 44)

Then, the third graph presents the discrepancy of other samplers. It can be seen that as expected, the Whitenoise sampler has a discrepancy in O(1/n), similarly to the discrepancy of regular and hexagonal grids (due to samples alignment). The N-Rooks sampler has an interesting behaviour, as its discrepancy seems to increase as the number of samples increase. That is most likely because N-Rooks shuffles a canonical arrangement to preserve its 1-D projections, and as the number of samples increases, samples are getting closer and closer to each other in those projections. Furthermore, in the 2-D plane, there is nothing that prevents samples from clumping. Thus, the discrepancy gets closer and closer from the discrepancy of the Whitenoise sampler as the number of samples increase. Another interesting point on this graph is the discrepancy of stratified sets, that is surprisingly good. Thanks to the stratification of samples, clumping is partly alleviated and aligning samples are less likely. Hence the improvements in the discrepancy. Finally, our last graph presents how our new samplers (described Chapter 5 on page 93 and Chapter 6 on page 107) behave regarding state of the art samplers. It can be seen that, as expected, they experimentally behave as Low Discrepancy Samplers. They thus present a better discrepancy than Blue Noise and Historical samplers.

However, Graph 7.5 on page 137 presents that, for our second method, the discrepancy is very degraded between the octaves. The size of such octaves depends on the subdivision factor K of our tile based system. A point set is therefore low discrepancy only if its size n verifies n = a 2 K , a ∈ N. Therefore, the higher the K, the more difficult to obtain a set that is Low discrepancy. On the other hand, the higher the K, the more Blue Noise our sets can be. This will be discussed further in the Chapter 8 on page 165.

l 2 star discrepancy

As the l ∞ star discrepancy is not tractable for large sets in more than 2-D, we need another measure to compare how the discrepancy of a sampler evolves as the dimension increases. To do so, we used the l 2 star discrepancy measure, as described Section 2.1.2 on page 19. This discrepancy was computed from the Generalized L2 Discrepancy tool from UTK. This tool takes as input a pointset P n , and computes its generalized l 2 star discrepancy using the following formula (previously given Equation 2.13)

(D 2 (P n )) 2 := 4 3 s - 2 n n i=1 s j=1 3 -x (i) j 2 2 + 1 n 2 n i=1 n i =1 s j=1 2 -max x (i) j , x (i ) j . (7.1)
This lets us compute this discrepancy with a O(sn 2 ) algorithm.

Similarly to the l ∞ star discrepancy, we first present the evolution of the l 2 star discrepancy for state of the art Low discrepancy samplers (Fig. 7.6 on the following page), then for Blue Noise samplers (Fig. 7.7 on the next page), and finally for other samplers (Fig. 7.8 on page 139).

Then, the discrepancies of the main samplers will be compared to the discrepancy of our new methods (Fig. 7.9 on page 139).

On all those graphs, we presented, for comparison, the theoretical convergence speeds of the l ∞ star discrepancy. Note that we do not expect the convergence speeds of the l 2 star discrepancies of samplers to follow those curves. However, it can still be noted that the measurements are quite close from those curves.

Graph 7.6 on the following page presents the l 2 star discrepancy of low discrepancy samplers. It shows that the convergence speed is still the same for all those samplers, and we see the same periodical increase in the discrepancy of Sobol [START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF] and Faure [START_REF] Faure | Discrépance de suites associées à un système de numération (en dimension s)[END_REF] in between the octaves, as was observed for the l ∞ star discrepancy.

For Blue Noise samplers, presented Graph 7.7, we note that we no longer observe the weird behaviours that we could see with the l ∞ star discrepancy. This is most likely due to the fact that the l ∞ star discrepancy keeps the maximal difference between box area and fraction of samples in this box (Section 2.1.1 on page 17). Thus, this discrepancy is highly sensitive to outliers, if a single measure fails, the whole discrepancy behaves abnormally. The l 2 star discrepancy on the other hand averages the differences for all boxes. Outliers are thus invisible in the final measure.

The last two graphs, presenting the l 2 star discrepancy of other samplers (Graph 7.8 on the next page), and comparing our methods to state of the art methods (Graph 7.9 on the facing page), lead to roughly the same conclusions as the ones from l ∞ star discrepancy graphs. Contrary to the analysis of 2-D point sets using the l 2 star discrepancy, for higher dimensional sets, we cannot present measurements over Blue Noise samplers as those samplers are of very limited dimensionality. Figure 7.10 thus only present the evolution of the l 2 star discrepancy of state of the art 4-D samplers compared with the 4-D discrepancy of our higher-dimensional method (Chapter 6).

We see in Figure 7.10 that, as expected, our projection scrambling does not preserve the full 4-D discrepancy. However, its discrepancy is still better than the discrepancy of a stratified sampler for lower number of samples. It seems that each stage of scrambling messes up the 4-D discrepancy of the set a little bit more, and that as the number of samples increases, our set has a discrepancy with a similar behaviour as the stratified sampler.

Integration error

In this section, we compare how various samplers perform when used for the Monte Carlo estimation of a given integrand.

On the following graphs however, we did not plot the integration error as

Δ = |I n -I|. (7.2)
Plotting this function for a single integrand would have led to very noisy curves and would have prevented an easy reading of the graphs. Instead, we used the theorem from [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF] that states that for unbiased stochastic samplers, their integration error is equivalent to the square root of the variance in integration error. Plotting this measure leads to much smoother curves, even for a single integrand. For deterministic samplers, we use the Cranley-Patterson shift (Section 4.1 on page 74) to randomize them, and also plotted their variance. This is possible since this Note however that integrating a single function with a Cranley-Patterson shifting of samples can be seen as integrating a Cranley-Patterson shifted function with a single set of samples. This may, in particular cases, create discontinuities in the function that were not present in the original integrand (Fig. 7.12). Graph 7.13 on the following page presents a comparisons of the MSE measured when integrating a cornered Gaussian function from a single Sobol and a single Stratified set and the MSE measured from averaging several Cranley-Patterson shifted realisation of those sets. Such integrand when shifted creates discontinuities which changes the frequency radial profile of the function. In all those integrands used for our following test, discontinuities were already present in the original integrand and thus, shifting it has no impact on the MSE. We observed that integration errors for functions with higher exponents (in absolute value) in their mean spectral profile are similar for all samplers. This is consistent with the results from [Pil+15], as a higher exponent means that the radial power spectrum from the integrand will be zero almost immediately and thus only the very low frequencies will be present. As a vast amount of samplers try to be as close to 0 as possible in the low frequencies the Fourier spectrum of all those samplers will have a similar impact on the observed error. Therefore, measuring a difference between those samplers for such integrands is difficult. The differences appear when integrating over functions having a mean spectral profile with a higher exponent.

We start by presenting this measured integration error when integrating an analytical 2-D disk and analytical 4-D spheres. Then, we present how samplers perform when integrating over a discretization of a natural scene, given as a HDR image. Finally, we present a new tool (still at an early development stage) to allow analytical integration of complex arrangements, that we hope may, with further work, allow to simulate the spectral properties of natural scenes, but with an analytical representation.

Integration over analytical functions

We start by integrating an analytical 2-D disk defined with

f Disk (x) := ⎧ ⎨ ⎩ 4 √ π if |x| < 0.25, 0 otherwise. (7.3)
This Disk is an interesting analytical integrand as it presents discontinuities which are very common in rendered scenes. One could also use a Triangle or a Square, but their Power spectra converge faster than the Power spectrum of the Disk, which converges in O(ρ -3 ). We can see on Figure 7.14 that not all Blue Noise samplers are equal for Monte Carlo estimators.

Consistently with what was observed for the discrepancy of the CVT [START_REF] Du | Centroidal Voronoi tessellations: Applications and algorithms[END_REF] and CapCVT [START_REF] Chen | Variational blue noise sampling[END_REF] samplers, their integration error suddenly increases drastically as the number of samples gets too high. Similarly for the discrepancy, we do not know at this time why this behaviour occurs but we assume it is either due to a misprogramming/misusage of the code (we are using the code from the original authors), or to a very minor shifting of the samples that may not be visible to the eye but that causes the integration error (and discrepancy) to increase drastically. It can also be noted that the Poisson Disk sampler, consistently with the results of [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF], performs quite poorly. Similarly, the FPO sampler has a higher integration variance. This is still unexplained but our hypothesis is that this is due to similar reasons as for the Poisson Disk sampler.

When looking at the integration error of Low discrepancy samplers, Figure 7.15 on the following page, we can note that all samplers performs very well. This is consistent with the Koksma-Hlawka theorem which bounds the convergence of the integration error of 2-D Low discrepancy samplers to O(log(n)/n).

The graph Figure 7.16 presents the results from other samplers. The behaviour of Whitenoise and stratified samplers are consistent with their theoretical convergence. However, the behaviours of NRooks and of the Hexagonal grid are surprising. For NRooks, this is most likely due to the regular canonical input that, once shuffled, does not allow for a covering of the final set as uniform as a White noise set. For the Hexagonal grid, the high variance observed within the variance values themselves prevents a proper analysis of what happens, but it seems to suffer from a different issue than the NRooks sampler as its discrepancy behaved perfectly well.

Finally, Figure 7.17 on page 145 compare the behaviour of our methods with respect to the main state of the art samplers. It can be seen that our integration error converges as fast as for usual Low discrepancy samplers, which was to be expected as our samplers are Low discrepancy. We also performed measurements for the integration error of a 4-D sphere, with 4-D samplers.

The results are given Figure 7.18 on the following page. It can be seen that, despite not being a Low Discrepancy 4-D sampler, our sampler performs very well, as, in the worst case, our set resumes as a 4-D stratified set which has a good performance regarding the integration error.

Integration over discrete natural scenes

In this section, we show our results when integrating over HDR images. We chose two HDR images with very different Fourier spectra profiles. Natural images spectra can be roughly categorized into two categories [START_REF] Torralba | Statistics of natural image categories[END_REF]: The ones showing man made scenes (buildings, streets, etc.) and the ones showing natural scenes (forests, grass, etc.). Man made structures have a much more anisotropic spectra and present a mean spectral profile in O(ρ -1.8 ) while natural images are much more isotropic with a mean spectral profile in O(ρ -2 ). The two gray scale HDR images are taken from the sIBL archive with a size of 1600x1600. Figure 7.19 on the next page presents those images along with their Fourier spectrum, radial profiles, and exponent.

The results from this section are similar to the results observed for analytical integrands. They are presented in Figures: 7.21, 7.20, 7.22, 7.23, 7.25 on page 149, 7.24, 7.26 and 7.27 on the following pages. Note that as the exponent of their radial profile is lower, distinction between samplers is made easier.

Integration over analytical arrangements

In this section, we present a new integration test. When integrating from natural HDR images, we only have so many images we can integrate from. Furthermore, those images are already discretized and we do not have any analytical expression for the integrand. Here, our final goal is to generate random synthetic analytical images, with a controllable profile spectrum, with arrangements of boxes and triangles.

We therefore generate randomly 2-D boxes (or triangles) of various sizes within a given domain.

By controlling the number of shapes of a given size, we can control the amount of higher/lower frequency in the final image. This control could be made finer by using white,blue,pink or green distributions of shapes to increase a given range of frequencies. As each shape is analytically integrable, we have the analytical integrand for the first shape, which is simply updated each time we add a new shape within the arrangement. Technically, we do this using the CGAL library [Cga], which provides efficient tools to handle such 2-D arrangements.

However, in our current implementation, to create integrands with a radial profile similar to natural images (between O(ρ -1. 5 ) and O(ρ -2 )), we need to generate very complex arrangements (containing around 15000 elements of various sizes and orientations (Fig. 7.28)). Unfortunately, when used for integration, such arrangements lead to biased results. We do not know for sure today what causes this bias but it seems that the analytical evaluation of those arrangements leads to numerical errors.

In this section we present integration error as measured for simpler and unbiased analytical arrangements. Figure 7.29 on the following page, with radial profiles between O(ρ -2 ) and O(ρ -2.5 ) (Fig. 7.29 on the next page).

The Box arrangement presents more anisotropy than the Triangle arrangement and emphasizes the behaviour of samplers when integrating axis aligned discontinuities. When using those synthetic images to measure convergence speed in integration, we get consistent results as those that would be obtained from true HDR images (Figures 7.30 on the facing page and 7.31 on the next page). Note that, as computing those arrangements takes some time, in this section we plotted the true MSE instead of the variance over several randomized realisations of the samplers/the integrands, hence the noisy curves.

Aliasing

In the previous sections, we studied how approximating a single integrand converges. However, this is insufficient in a rendering context, where many integrands are solved (one for each pixel), and where the set of all approximations values are then reconstructed into a full image. We study in this section how the integration error from each samplers is distributed within this full image. A stochastically distributed error leads to noise, and a regularly distributed error leads to aliasing.

We present here results from reconstructing a 2 -D function with a wide range of frequencies, and results of various samplers in a true rendering context.

Zoneplate

Our first test to measure this is the reconstruction of the 2-D Zoneplate function. This reconstruction is done by sampling the function with a 2-D sampling pattern and then using a Mitchell-Netravali filter [START_REF] Don | Reconstruction filters in computer-graphics[END_REF] to retrieve the original function from all the measured discrete values. See Appendix A on page 168 for more informations. This Zoneplate function f Z is illustrated Figure 7.32 and is formally expressed by f Z (x, y) := sin(x 2 + y 2 ). (7.4) It is interesting to us as this function presents a very wide range of frequencies, which is very difficult to capture when discretizing. Aliasing becomes obvious when using this function, and its relationship to the Fourier spectrum is also easily visible. The zoneplate images from Figure 7.33 show that Blue Noise samplers allow for a good reconstruction up to a certain frequency. Above that frequency, the reconstruction turns to noise. Among this noise, an arc of stronger noise artefacts appears. This arc is due to the peak in the Blue Noise fourier spectrum. Similarly, in Figure 7.34, the Low discrepancy samplers would perform a very good reconstruction, as no noise is visible. The issue is that instead, it turned this noise into very visible artefacts. Those artefacts actually align with the peaks in the Fourier spectrum of the sampling pattern used to do this reconstruction. Finally, the samplers presented in Figure 7.35 usually generate more noise than Blue Noise sampling patterns, and the lattice samplers present artefacts stronger than the ones from the Low discrepancy samplers.

In Figure 7.33, the zoneplate reconstruction from the samplers of our methods are presented. Those reconstructions are very similar to the ones from classical Blue Noise sampler, which was to be expected as our methods present a Blue Noise Fourier spectrum.

Renderings

In this section, we present the performances of our second method (Chapter 6 on page 107) in a rendering context. All our rendering examples were performed using the PBRT engine [START_REF] Pharr | Physically based rendering: From theory to implementation[END_REF]. For this whole section, we strongly advise the reader to refer to the pdf version of this document as the printing may have degraded the quality and readability of the results.

To use our sampler within this engine, we generated a set of N s-D samples (usually with N = 2 18 or N = 2 20 ), and tiled this set so that each pixel is covered by the expected number of samples. We could have instead generated independently but this would have led to correlation

Aliasing

between pixels that would have led to very disturbing artefacts. We chose to tile our sampler as, for implementation simplicity, we are reading the samples from an input file instead of having them generated on the fly by the system. Thus, to avoid reading and storing huge files, we chose to tile a smaller set instead. We also used this method to perform renderings from Owen's scrambled Sobol sets.

We compared our sampler to Sobol and Halton, which are both generated using the same global stretching of N samples. The only difference being that as Sobol and Halton are defined analytically, it is possible to generate billions and billions of samples and thus no tiling is required. On the contrary, the Stratified sampler renders a scene by generating a new independent set of samples within each pixel.

We rendered several scenes, and compared the mean squared error (MSE) of the rendered scene to the ground truth scene. Note that the MSE given here is not the same as the MSE given in Section 7.2 on page 140. Section 7.2 on page 140 measures the error/variance when integrating a single scene, returning a single value as the integrand of the whole scene. The rendered scenes here are made from an aggregation of approximated integrand values, one for each pixel. The MSE measured here can therefore be perceived a sum of all the MSE from all the integration errors from all pixels. It does not take into account how the error is distributed, and is not affected by aliasing artefacts.

In Figure 7.37, we start by using our sampler solely for 4 dimensions. This scene uses 8 dimensional samples: 2 for the image plane, 2 for the point on the light, 2 for the light scattering, and 2 for the object reflectance function. We rendered it using various samplers to sample the 2D image plane and the 2D light source, while all other dimensions sampled using the Sobol sequence. Since the Sobol sequence is deterministic, differences between the images are only due to the difference in sampling for the image and light and therefore allow for an easy comparison of the different samplers. It can be seen that our sampler has a MSE similar to Sobol, but no longer presents aliasing in the shadow. Note that we used our sampler with 64 samples per tile (K = 8), and the Sobol indices (1, 4) and (2, 3).

In Figure 7.38, we render a complex scene with various samplers. In the case of our sampler, we used 64 (K = 8) samples per tile and optimized 3 2-D projections, suing the Sobol indices (1,4), (2,3) and (5,7). Those optimized projections were used for the image plane, the light and the lens (for the depth of field effect). We used our hierarchical Owen' scrambling (Section 4.2 on page 79 for the remaining dimensions. Similarly to Figure 7.37, we note that our approach's MSE is better than most methods as is is similar to Sobol's. However, we present much less aliasing.

Figure 7.39 presents the difference in renderings when using 16 (K = 4) or 64 (K = 8) samples per tile on a very complex scene requiring 20-D samples. We used the same setup as for Figure 7.38, which is optimizing the projections for the image plane, the lens, and the light and filling the other dimensions with Owen' scrambling. Note that when using 4 samples per pixel, the crown appears darker. This phenomenon is still unexplained.

Timings

In this section we present timings measured from the samplers implementation from the UTK framework. Table 7.1 presents the times in seconds need for a sampler to generate a given amount of samples. it can be noted that the analytical Low discrepancy samplers are very efficient, when the Blue Noise samplers based on optimization cannot generate more that 10 5 samples. All those tests were performed on a CPU Intel® Core™ i5-4440 CPU @ 3.10GHz × 4.

Conclusion

The experiments presented in this chapter led to the following conclusions. First, they experimentally confirm our intuition regarding the Koksma Hlawka theorem. Meaning that despite that our scenes are not of finite variations in the sense of Hardy and Krause (as required by the original Koksma Hlawka theorem), the discrepancy seem to still be related to the integration error as an upper bound. We observe thar the Low discrepancy samplers present an integration error that is related to their discrepancy, and decreases in O(log(n)/n) (in 2-D). It also confirmed the theoretical convergence speeds of various stochastic samplers from [START_REF] Pilleboue | Variance Analysis for Monte Carlo Integration[END_REF] and listed in Table 3.1 on page 70. The Zoneplate test allowed to see the strong aliasing artefacts that occurs when peaks are visible in the Fourier spectrum, and the renderings allow to see how a structured sampler leads to a structured repartition of the error in rendering. For all those experiments, our samplers performed as expected. Meaning that their Low Discrepancy gives them a convergence speed in O(log(n)/n) (in 2-D), and their Blue Noise spectrum removes all structured error from the renderings. This confirms that our new proposed solutions, contrary to previous sampling patterns, performs well both for aliasing and integration error.

Conclusion

Conclusion

Generating a photo-realistic image on a computer is a complex problem, that involves solving high dimensional recursive integrands. Unfortunately, those integrands cannot be solved analytically. Instead, we rely on a numerical approximation of the integrand, evaluated using a Monte Carlo estimator. Roughly, such estimators sample the integrand, generating several discrete values, and deduce the integration result as a sum of all those values. There are several ways to chose which sample of the integrand should be taken, each leading to different results. In this PhD, we have studied many sampling patterns and evaluated how they perform in a rendering context. A sampling pattern is considered a good sampling pattern if the Monte Carlo approximation converges as fast as possible to the true integration value, and if the integration error is not structured in the final image (if there aren't any aliasing artefacts).

It was shown that samplers that are said to be Low Discrepancy are the ones with the fastest convergence (in the worst theoretical case). However, they usually also come with a very regular structure, which also structures the error in the final image. Furthermore, the Low discrepancy guarantees a uniform coverage of the domain only asymptotically. As example, the Sobol sequence is Low discrepancy, but may present a distribution of samples aligned on the diagonal when the number of samples taken is too low.

On the other hand, the better non-aliased sampling patterns are the Blue Noise patterns. Such samplers ensure a better covering of the domain by forcing all samples to be equidistant from each other. However, they asymptotically converge slower than Low Discrepancy samplers.

Our contributions aimed at developing samplers that would be Low Discrepancy, and present a Blue Noise spectrum. We demonstrated that from this association we have sampling patterns that have an optimal convergence speed in integration, and an absence of aliasing in the final rendered image. Furthermore, compared to pure randomization of a low discrepancy set, the Blue Noise spectrum improves the approximation by ensuring a more uniform covering of the domain, even for a given low number n of samples, where the asymptotic low discrepancy property cannot apply.

Our first contribution (Chapter 5 on page 93), aimed at providing the theoretical and empirical proof that a Blue Noise spectrum is not incompatible from a Low Discrepancy property. However, this contribution is limited to 2-D point sets.

Our second contribution (Chapter 6 on page 107), allows to scramble specific 2D projections of a Sobol sequence to generate higher dimensional patterns with Blue Noise and Low discrepancy properties in those projections. Therefore, during rendering, the integration error decreases quickly, and no aliasing appears in the final image. Furthermore, our method generates sequences of samples.

However, this contribution still leaves the path open future work. First, the K factor is a difficult trade off. As noted in Section 7.1.1 on page 133, our projections are only low discrepancy for sets of size n = a 2 K , a ∈ N. Therefore, the higher the K value, the more degraded the discrepancy. But on the other hand, with a higher K, we can hope for a better spectrum for the final set. An aware reader may have already noted that this drawback is due to the fact that our method generates (0, a, 2)-nets in base 2 K . This choice of a higher basis was motivated by the intuition that that the Sobol sequence is the only existing (0, 2)-sequence in base 2 (modulo the order of the samples). As I could not find formal proof of this, this is still conjecture. Nonetheless, creating Blue Noise sequences from digital nets seems to require looking among sequences defined in basis b > 2.

Still, the main reason we aim at having a Blue Noise spectrum comes from the fact that, for a long time, aliasing was considered worse that stochastic noise as it is much more noticeable to the human eye. However, such artefacts only occur for non-converged renderings, and such renderings are most the time used as input for denoising algorithms (such as the one from [START_REF] Schied | Spatiotemporal variance-guided filtering: real-time reconstruction for path-traced global illumination[END_REF]). Such algorithms denoise images generated with as little as 1 sample per pixel. Therefore, choosing the type of noise to generate in an image should maybe be weighted against what denoising algorithms can efficiently remove, instead of what the human eye perceives.

With their structured error, aliased low discrepancy samplers may be in fact easier to denoise than stochastic samplers.

Finally, in the particular case of rendering, many approaches rely on importance sampling to adapt the density of samples pixel-wise, depending on the complexity of the scene behind those pixels. Such methods have been further improved through the use of deep learning, leading to impressive results [START_REF] Dahm | Learning Light Transport the Reinforced Way[END_REF]. In our contributions, we either have no adaptivity, or an adaptivity constrained by the precomputed optimized sets within each tile. To accommodate importance sampling, such samplers need to be made more adaptive, maybe using optimization based techniques as the one presented in Section 4.3 on page 80, assuming such methods could be performed in real-time.

Tutu

A Discretization and Reconstruction of a continuous signal

In this section, we present the theoretical basis to reconstruct a signal from a set of evaluations. We will start by defining what is a signal, and how we can discretize a continuous signal. Then, we will see under which conditions we can properly reconstruct the original signal from its discrete form.

All the notations we use regarding signal theory follow the ones from [Gla95].

⊗ ω = ω F f F III T = III 1 T F f ⊗

A.2 Spatial and Frequency spaces

Any CT or DT signal can be expressed with many different forms. The ones that will interest us here are its spatial-representation and its frequency-representation. Those forms are strictly equivalent, and represent the same signal in a different space (space is to be understood here as a choice of referential basis and not as a geometric space). The spatial-representation of a signal is its expression in the spatial-space, which is the most common one. The frequency-representation of a signal is its expression in the frequency space, where the referential is composed of complex exponentials [START_REF] Newbold | The Fourier transform and its applications[END_REF].

The frequency representation of a signal is often referred to as its Fourier spectrum. To transform a signal from its spatial to its frequency representation, we use the Fourier transform. We denote F f , the frequency representation of a spatial signal f , obtained from computing the Fourier transform of f . The inverse operation can be done using the Inverse Fourier transform. We denote F -1 g the spatial representation of a frequency signal g, obtained from computing the Inverse Fourier transform of g. Note that f = F -1 F f . We will not get into much more detail here as this is explained thoroughly Section 2.2.1 on page 20.

Even though those representations are strictly equivalent, applying an operator onto a signal in spatial space is not equivalent to applying the same operator to this signal in frequency space.

The property that interests us here is that a product applied in the spatial space is equivalent to applying a convolution in the frequency space, and vice versa. More formally,

f • g = F f ⊗ F g F f • F g = f ⊗ g, (A.5)
where ⊗ denotes the convolution operator. This convolution operator performs the product of every value of a function f with every value of a function g

ω × ω = ω F f F sinc = B H F f • B H
(f ⊗ g) (t) := Θ f (θ)g(t -θ)dθ, (A.6)
where Θ is the support for the function f . Furthermore, we will note that we have the following properties for the Dirac impulse and Dirac comb: a Dirac impulse in spatial space transforms to a Dirac impulse in frequency space, and similarly, transforming a Dirac comb only inverts the frequency of its set of impulses. More formally,

δ = F δ (A.7) III T = F III 1 T . (A.8)
This means that computing the product of a CT signal f with a Dirac comb III T in the spatial space (to discretize it) resumes in the frequency space as a convolution of the Fourier transform of f with the Fourier transform of III T . One of the consequence of this operation is that it repeats the signal creating aliases of the original signal in frequency space (Fig. A.5 on the preceding page). It can be seen that in frequency space, if the frequency of the Dirac comb is too high (which means, if the frequency of the Dirac comb is too low in spatial space), the aliases will cover each other, leading to a different representation of the signal in the frequency space (Fig. A.10 on page 173). When reconstructing our signal, this will have a huge impact.

A.3 Reconstruction

In spatial space, we can reconstruct our signal by convolving it with a chosen reconstruction filter. This means, computing the product of the frequency representation of our signal with the frequency representation of our filter. As computing this product is much easier than performing a convolution, we will focus on what happens in the frequency space.

We recall that we have a discretized version of our signal in frequency space, which contains several aliases of the frequency representation of the original signal (Fig. A.5 on the preceding page). Our goal is therefore to fully retrieve the original alias and only the original alias. The easiest way to do so is to compute the product of our signal in frequency space with a box function B H defined as: In this appendix, we present a contribution we made in a very different field. During this PhD, we focused mainly on solutions to numerically approximate the rendering equation. However, we also worked on solutions to analytically simplify this equation, in particular cases when the surface follow a given statistical distribution. In those cases, the integration sums up as the CDF of the distribution characterizing the surface. Such methods did not give the expected results, but we were able to reuse this work for the particular case of digital curvature estimation.

t ⊗ t = t f sinc f (t) ⊗ sinc
B H (t) = ⎧ ⎨ ⎩ 1 H , if t ∈ [-H 2 , H 2 ] 0, otherwise .
In this appendix, we introduce what is curvature estimation of digital objects, and how thanks to the capacity of modern GPUs we could perform this computation in real time. Along with this, we also perform an on the fly triangulation of the digital object, with proper level of details, to ensure a good rasterization of the triangles (where there never are two triangles projecting into the same pixel).

B.2 Introduction

Volumetric objects are being more and more popular in many applications ranging from object modeling and rendering in Computer Graphics to geometry processing in Medical Imaging or Material Sciences. When considering large volumetric data, interactive visualization of those objects (or isosurfaces) is a complex problem. Such issues become even more difficult when dynamic volumetric datasets are considered. Beside visualization, we are also interested in performing geometry processing on the digital object and to explore different parameter settings of the geometry processing tool. Here, we focus on curvature tensor estimation (mean/Gaussian curvature, principal curvatures directions. . . ). Most curvature estimators require a parameter fixing the scale at which the computation is performed. For short, such parameter (integration radius, convolution kernel size. . . ) specifies a scale for analyzing the object surface, and is naturally related to the amount allowed perturbations on input data. As a consequence, when using such estimators, exploring different values of this parameter is mandatory. However, this is usually an offline process, due to the amount of computations that need to be done.

Contributions

In this work, we propose a framework to perform interactive visualization of complex 3D digital structures combined with a dynamic curvature tensor estimation. We define a fully data parallel process on the GPU (Graphics Processor Unit) to both efficiently extract adaptive isosurface and compute per vertex curvature tensor using Integral Invariants estimators. This system allows us to visualize curvature tensor in real-time on large dynamic objects. Our approach combines a GPU implementation of pointerless octrees to represent the data, an adaptive viewpoint-dependent mesh extraction, and a GPU implementation of integral invariant curvature estimators.

Related works

Extracting and visualizing isosurface on volumetric data has been widely investigated since the seminal Marching Cubes approach by Lorensen and Cline [START_REF] William | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF]. This approach being data parallel, GPU implementation of this method is very efficient [START_REF] Tatarchuk | Real-time isosurface extraction using the GPU programmable geometry pipeline[END_REF]. However, such technique generates a lot of triangles which is not well suited to large digital data. Hence, adaptive approaches have been proposed in order to optimize the triangulated mesh resolution according to the object geometry or the viewpoint. Curvature estimation on discrete or digital surface has also been widely investigated.

In [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF], authors propose digital versions of Integral Invariant estimators [START_REF] Pottmann | Principal curvatures from the integral invariant viewpoint[END_REF][START_REF] Pottmann | Integral invariants for robust geometry processing[END_REF] in order to estimate the complete curvature tensor (mean/Gaussian curvatures, principal curvatures, principal curvature directions, normal vector field) on digital surfaces. Such approaches are based on an integration principle using a ball kernel of a given radius. Additionally, authors have demonstrated that these estimators have multigrid convergence properties. In this article, we also propose an efficient GPU implementation of such estimators to visualize such curvature fields in real-time and to interactively change the value of the kernel radius.

In this paper, we first present the previous works on curvature tensor estimation. Then, we propose an efficient GPU approach to extract a triangulated isosurface from a digital object. Finally, we propose a fully-parallel GPU pipeline to compute curvature tensor in real-time and present our results.

B.3 Curvature Tensor Estimation

In our context, we consider digital shapes (any subset Z of Z d ) and boundaries of digital shapes Bd(Z). We denote by G h (X) the Gauss digitization of a shape X ⊂ R d in a d-dimensional grid with grid step h, i.e. In [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF], authors define a 2D digital curvature estimator κR and a 3D digital mean curvature estimator ĤR based on the digital volume estimator Vol(Y, h) := h d Card(Y ) (area estimator Area(Y, h) in dimension 2):

Definition 1 Given the Gauss digitization Z := G h (X) of a shape X ⊂ R 2 (or R 3 for the 3D mean curvature estimator), digital curvature estimators are defined for any point x ∈ R 2 (or R 3 ) as: for convex shapes with at least C 3 -boundary and bounded curvature (setting R := kh 1 3 for some k ∈ R).

∀0 < h < R,
In [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF], authors have also defined 3D digital principal curvature estimators κR 1 and κR 2 on Z ⊂ Z 3 based on digital moments: Definition 2 Given the Gauss digitization Z := G h (X) of a shape X ⊂ R 3 , 3D digital principal curvature estimators are defined for any point x ∈ R 3 as: where λ1 and λ2 are the two greatest eigenvalues of the covariance matrix of B R/h (x/h) ∩ Z.

∀0 < h < R,
The covariance matrix needs to compute digital moments of order 0, 1 and 2 (see Equation 20 of [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF] for more details). These estimators are proven convergent in O h , where k is a constant related to the maximal curvature of the shape, on convex shapes with at least C 3 -boundary and bounded curvature [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF]. Additionally, eigenvectors associated to λ1 and λ2 of the covariance matrix are principal curvature direction estimators ŵR 1 and ŵR 2 . The smallest eigenvector corresponds to the normal direction nR at x. Convergence results can be found in [START_REF] Lachaud | Robust and Convergent Curvature and Normal Estimators with Digital Integral Invariants[END_REF].

It has been shown that the radius of the ball depends on the geometry of the underlying shape. In [START_REF] Levallois | Parameter-free and Multigrid Convergent Digital Curvature Estimators[END_REF], a proposal was made for a parameter-free estimation of the radius of the ball by analyzing the shape w.r.t. the local shape geometry using maximal digital straight segments of the digital boundary. In [START_REF] Levallois | Scale-space Feature Extraction on Digital Surfaces[END_REF], these estimators have been analyzed in scale-space (for a range of radii) for a given digital shape. This allows to detect features of the shape thanks to the behavior of estimators on singularities. As a consequence, for all these integral invariant based approaches, we need to consider different ball radius which could be time consuming when implemented on CPU. We propose here a fully parallel implementation on GPU allowing us to change the radius and thus update the estimated quantities in real-time.

B.4 Isosurface Extraction on GPU

In this section, we detail the adaptive isosurface extraction algorithm. The proposed approach uses an octree representation of the input object on which an adaptive Marching Cube builds the isosurface efficiently. Such hierarchical representation of the object allows us to handle large datasets and to locally adapt the level of details w.r.t. the geometry or camera position. We first present the octree representation and then the isosurface extraction. Representing a hierarchical structure on GPU is usually challenging since such a data parallel component is unable to handle recursivity. Efficient spatial tree encoding can be achieved using pointerless structures such as linear quadtrees or octrees Gargantini A cell's code encodes its position w.r.t. its parent cell and its complete path to the tree root. Hence, the tree is fully represented as a linear vector of its leaves. Furthermore, cell operations such as subdivision, merging can be efficiently implemented using bitwise operations on the Morton code. In the following, we use the GPU friendly implementation proposed by Dupuy et al. [START_REF] Dupuy | GPU Pro 5[END_REF].

B.4.1 Linear Octree Representation

B.4.2 Data parallel and adaptive mesh generation

Using this spatial data structure, a triangulated mesh can be constructed using Marching Cubes [START_REF] William | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF] (MC for short): the triangulation is generated from local triangle patches computed on local cell configurations. Such approach is fully parallel and easy to implement on the GPU. However, since adjacent cells may not have the same depth in the octree, original Lorensen and Cline's rules need to be updated (see In the following, we use the algorithm proposed by Lengyel et al. [START_REF] Stephen | Voxel-based terrain for real-time virtual simulations[END_REF]. First, this approach forces the octree structure to make sure that the depth difference between any two adjacent cells is at most one. Then, Lengyel et al. introduce the concept of transition cells. Those cells are defined to be inserted between two neighboring octree cells of different depth. With specific MC configurations for their triangulation, a crack free mesh can be extracted.

Similarly to original MC algorithm, this approach is well suited to a GPU implementation: given a set of cells (a vector of morton codes), each triangle patch can be extracted in parallel for both regular cells and transition cells.

B.4.3 Level of Details Criteria and Temporal Updates

As illustrated in Figure B.2 on the previous page-middle, we propose a viewpoint dependent criterion to decide if a cell needs to be refined: the closer we are to the camera, the finer the cells are. Such a criterion is also well suited to GPU since it can be evaluated independently on every cell. Our distance criterion is based on the ratio (visibility ratio in the following) between the cell diameter l(c) (power of 2 depending on the depth), and its projected size. For a given cell c, split and merge decision are based on this visibility ratio:

• c is split if its children cells have a visibility ratio greater than constant k ;

• c and its sibling cells are merged if their parent cell c has a visibility ratio lower than k ;

• otherwise, the cell c stays for the next frame.

Using such a criterion, split and merge decisions are computed in parallel from the morton codes of all cells. Once decisions have been made, a new set of cells is sent to the mesh generation step described above. Finally, before constructing the triangulation from remaining cells and transition cells, geometrical culling is performed in order to skip the triangle patch construction for cells that are not visible.

B.5 Interactive Curvature Computation on GPU

We first present design principles for the GPU implementation and then present our approaches. The general idea is to perform an Integral Invariant computation on GPU at each vertex of the generated triangulated mesh. Since GPU have massively parallel architectures, we can do all those computations in parallel to obtain a very efficient implementation. Please note that when the LoD criterion is removed, each MC vertex is exactly centered at a surfel center. At different depth of the octree, MC vertices still correspond to surfel centers of subsampled versions of the input object. Hence, Integral Invariant framework defined in Section B. In this case, errors only occurs for cells lying inside X (with density 1) not entirely covered by B R (x). Furthermore, the mipmap texture can be used to design, using a single texture probe, a fast inclusion test of a given cell c at level l into the shapes: we say that c is in X if its density (retrieved from the texture probe at level l) is greater than 1/2. The idea here is to mimic a kind of adaptive Gauss digitization process.

Hierarchical Decomposition. Using the hierarchical nature of the mipmap texture, we could also consider hierarchical decompositions of B R (x). The idea is to decompose the ball into mipmap cells of different resolution in order to limit the number of texture access (important bottleneck on GPU hardwares) and to get better approximated quantities. On the GPU, computing a hierarchical decomposition of B R (x) at each point x is highly inefficient since the hardware optimizes the parallelism only when the micro-program described in the shader has a predictive execution flow. Hence, implementing the recursive algorithm (or its de-recursified version) requires a lot of branches (conditional if tests) in the flow. We have thus considered a fast multi-resolution approximation as illustrated in Fig. B.5-middle: for a given radius R, we precompute a hierarchical decomposition of B R . Such decomposition is made of octree cells at different resolutions. At a given MC vertex x, the algorithm becomes simple since we just scan the ball cells and estimate the area (or other moments) from the mipmap values associated to this cell. Note that these precomputed cells in the B R decomposition may not be aligned with mipmap cells. However, we can use GPU which interpolates mipmap values if we probe at non-discrete positions. For a given radius R, the expected number of cells in B R is in O(log R). Implementation details on the hierarchical octree representation can be found in the supplementary material (see Section B.7 on page 183).

Dynamic Refinement. In this last approach, we want to optimize the interactivity and the execution flow or parallelism on GPU. The integration is simply computed using a regular grid at different resolution l (Fig. B.5 on the facing page-right). The shader code becomes trivial (simple triple loop) and a lot of texture fetches are involved, but interactivity can be easily obtained. Indeed, we consider the following multi-pass approach:

1. When the surface geometry has been computed, we set the current level l at an high level l := l max of the mipmap pyramid.

2. We compute the integrals (and the curvature tensor) at the mipmap level l and send the estimated quantities to the visualization stage.

3. If the user changes the camera position or the computation parameters, we return to step 1.

4. If the current framerate is above a given threshold, we decrease l and return to step 2 if the final l 0 level has not been reached yet.

A consequence of the multi-pass approach is that when there is no interaction (camera settings, parameters), the GPU automatically refines the estimation. Even if step 2 is quite expensive, in

O R 2 l

3

, the interactive control in step 3 considerably improves the interactive exploration of the tensor with fast preliminary approximations which are quickly refined.

Compared to the hierarchical approach, no precomputation is required for a given radius R.

As a consequence, we could even locally adapt the ball radius to the geometry (for instance following the octree depth of the current MC vertex). In the next section, we evaluate the performances of both approaches.

B.6 Experiments

B.6.1 Full resolution experiment

We first evaluate curvature estimations on a full-resolution geometry obtained by disabling the LoD criterion at the mesh generation step. sults of curvature tensor estimation (mean, Gaussian, first and second principal directions) on "OctaFlower" at level l 0 (considered as our ground truth in the following). Then, we compare the approximations made by computing curvature from level l ≥ l 0 in Figure B.7 on the facing page. We can see that results using approximation seems to quickly converge to ground truth results. Table B.1 on the next page shows numerical results of L ∞ error (maximal absolute difference for all vertices) and L 2 error (mean squared errors of all vertices), as well the number of mipmap texture fetches. We can also see that the number of texel fetch, required to compute the curvature, reduces drastically when computing approximations. However, at higher levels, approximation errors become more important. Those levels are thus never used in practice, they are presented here to illustrate how the refining process converges to the l 0 level.

We also compare them with the hierarchical algorithm (as discussed in Section B.5.1 on page 180). This method introduces a higher error when compared with l 0 . This is due to the precomputation of the subdivision that no longer ensures to fetch data at the center of a mipmap cell. The GPU interpolation reduces the bias, but it does not remove it.

B.6.2 Fully adaptive evaluation

When dealing with large datasets, we cannot expect real-time curvature tensor estimation if we consider the full resolution geometry, due to the huge amount of data to process. By computing a dynamically refined approximate curvature tensor estimation joined with an adaptive triangulation (as discussed in Section B.5.1 on page 180), we manage to maintain a real-time framerate by giving control over the amount of data to process at each frame. In Figure B.8 on page 184, we compare timings (in logarithmic scale) for a triangulation that is dynamically refined according to its distance to the camera. We measured those timings using the multiresolution regular and the hierarchical algorithm. First we can note that the required time to compute the ground truth curvature for an object is usually as high as the time required to extract its geometry and to compute all of the above levels. Using approximations is thus mandatory. Another advantage with approximations is that it allows us to get a visualization of our object as soon as we run the application. This allows for real-time interactions with the object, required in order to change the visualized quantity, curvature radius, etc.

It is also visible in Figure B.8 on the next page that a hierarchical decomposition greatly reduces the curvature computation time, especially with big radii. However, due to the precomputation and the current hierarchical structure, this algorithm is biased and creates an error (presented in Table B.1) that needs to be considered. exploration in real-time on large datasets: "XYZ-Dragon" (with digital domain of 512 3 ) and "Snow microstructures" (233 3 ).

B.7 Conclusion and Discussion

In this article, we have proposed a fully data parallel framework on GPU hardware which combines an adaptive isosurface construction from digital data with a curvature tensor estimation at each vertex. Using this approach, we can explore in real-time different curvature measurements (mean, Gaussian, principal directions) with different ball radii on potentially large dynamic dataset. Our proposal relies on both a linear octree representation with Morton codes and an efficient integral computation on GPU. The source code and additional material (video, . . . ) are available on the project website (https://github.com/dcoeurjo/ICTV). List of Figures

FigureFig. 0 .

 0 Fig. 0.1.: (a) Présente l'ensemble de points 2-D initial de notre méthode. Cet ensemble est permuté sur chaque coordonnées (b,c) Pour donner un ensemble qui soit bruit bleu à l'arrivée (d). On peut noter que l'ensemble de départ a un spectre de Fourier très différent de celui de l'ensemble d'arrivée. On peut également noter que le spectre de Fourier de l'ensemble d'arrivée présente une croix en son centre. Cette croix est symptomatique de la basse discrépance et n'est en aucun cas un handicap aux capacités de notre échantillonneur en terme d'intégration.

Tout ceci est illustré Figure 0 Fig. 0. 2 .

 02 Fig. 0.2.: (a) Présente l'ensemble de points 2-D initial, répartis dans des tuiles à un niveau l. On permute individuellement chaque tuile (b,c) pour obtenir un ensemble qui soit bruit bleu à l'arrivée (d). Une tuile contient toujours le même nombre de points. Pour rajouter des points, on subdivise les tuiles. Cette subdivision peut être faite localement pour faire varier la densité de points (e).

(

  My Little Pony : Friendship is magic) I love you guys ! PS: Ah oui, je remercie aussi tous les chats qui m'ont ronronné dessus !

Fig. 1

 1 Fig. 1.1.: The rendering process.

  Fig. 1.2.: Different artefacts observed when rendering a scene using different sampling strategies.

  Fig. 1.3.: Various 2-D point sets P 1024 generated from different samplers S.

Fig. 2 .

 2 Fig. 2.1.: Spectral test in 2-D on Korobov lattices. Figure from [Lem09]

Fig. 2 . 2 .:

 22 Fig. 2.2.:Congruent boxes in 2-D on the Sobol sequence[START_REF] Il | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]

PFig. 2 . 3 .Fig. 2 . 4 .Fig. 2 . 5 .:Fig. 2 . 6 .:

 23242526 Fig. 2.3.: The CDF of a pointset P n

Fig. 2 .

 2 Fig. 2.8.: Fourier Serie decompositions F N of the square wave function and their corresponding Fourier transforms.

Fig. 2 .

 2 Fig. 2.9.: Power spectrum and Radial average characterization.

Fig. 2 .

 2 Fig. 2.10.: Illustration of Ripley's K function and of the PCF

Fig. 2 .

 2 Fig. 2.11.: Ripley's K-function, the L-function and the PCF of the pointset P n .

Fig. 2 .

 2 Fig. 2.12.: Product of the radial spectra of the sampling pattern Ps and of the integrand PI .

Fig. 2 .

 2 Fig. 2.13.: Blue Noise fourier spectrum and radial average.

  /(N sqrt(N)) Error^2 = k_2(log(N)/N)^2

Fig. 2 .

 2 Fig. 2.14.: Comparison of convergence speeds.

Fig. 3 .

 3 Fig. 3.1.: Characteristics of a Whitenoise sampler.

Fig. 3 . 2 .:

 32 Fig. 3.2.: Comparison between a Whitenoise point set and a Stratified point set.

Fig. 3 . 3 .:

 33 Fig. 3.3.: Characteristics of a Stratified sampler.

Fig. 3 . 4 .

 34 Fig. 3.4.: 2-D point set with stratified 1-D projections.

Fig. 3 . 5 .

 35 Fig. 3.5.: Shuffling canonical arrangements

Fig. 3 . 6 .:

 36 Fig. 3.6.: Characteristics of a N-Rooks sampler [Shi+91]

Fig. 3 . 7 .:

 37 Fig. 3.7.: Characteristics of a Correlated Multi Jitter sampler [Ken13].

Fig. 3 .Fig. 3 .

 33 Fig. 3.8.: Pointset and Radial Spectrum from several colors of noise.

Fig. 3 .

 3 Fig. 3.10.: Characteristics of a Relaxation Dart Throwing sampler [MF92].

Fig. 3 .

 3 Fig. 3.11.: Characteristics of Real Time Accurate Poisson Disk sampler [DH06].

Fig. 3 .

 3 Fig. 3.12.: Product of the radial spectrum of a Poisson sampling pattern Ps and of an integrand PI .

Fig. 3 .

 3 Fig. 3.13.: Characteristic of a Lloyd's relaxation [Du+99].

Fig. 3 .

 3 Fig. 3.14.: Lloyd's algorithm.

Fig. 3 .

 3 Fig. 3.15.: Characteristics of the Capacity Constraint Voronoi sampler [Bal+09]

Fig. 3 .

 3 Fig. 3.16.: Characteristics of the FPO sampler [Sch+11]

Fig. 3 .

 3 Fig. 3.17.: Characteristics of the BNOT sampler [DG+12].

Fig. 3 .Fig. 3 .

 33 Fig. 3.18.: Single iteration of a naive basic particle based Blue Noise sampling.

Fig. 3 .

 3 Fig. 3.20.: Wang Tiles paving.

Fig. 3 .Fig. 3 .

 33 Fig. 3.21.: Sampling patterns from [Ost+04],[START_REF] Ostromoukhov | Sampling with polyominoes[END_REF] and[START_REF] Wachtel | Fast tile-based adaptive sampling with user-specified Fourier spectra[END_REF] and their underlying paving system.

Fig. 3 .

 3 Fig. 3.23.: Characteristics of the AA Pattern sampler [Ahm+15].

Fig. 3 .Fig. 3 . 3 . 2

 3332 Fig. 3.24.: Halftoned images from various dithering algorithms.

Fig. 3 .

 3 Fig. 3.26.: (left) Stippling with isotropic Blue Noise [Bal+09] and (right) with anisotropic Blue Noise [Li+10a]. Figure courtesy of [Li+10a]

Fig. 3 .Fig. 3 .

 33 Fig. 3.27.: Characteristics of the Step Blue Noise sampler [Hec+13].

Fig. 3 .

 3 Fig. 3.29.: Korobov point set with a = 7 [Kor59]. Zoom to see the peaks in the Fourier spectrum.

Fig. 3 .

 3 Fig. 3.30.: (0, 4, 2)-net in base 2.

Fig. 3 .

 3 Fig. 3.31.: (3, 4, 2)-net in base 2.

3. 3

 3 Low discrepancy samplers Definition 4 ([Lem09, Def 5.7]) A (t, s)-sequence is a sequence of points for which each bary segment of the form v

Fig. 3 .

 3 Fig. 3.32.: Halton point set [Hal64].

  27) where I is the identity matrix, x = (.a 0 a 1 . . . a N ) b and i = (b N . . . b 1 b 0 ) b . The Van der corput sequence is often denoted φ.

Fig. 3 .

 3 Fig. 3.33.: Sobol point set with t = 0 [Sob67].

  x + 1, corresponding to indices 0 and 1. (a) Polynomials x 2 + x 1 + 1 and x 3 + x + 1, corresponding to indices 3 and 4.(a) Polynomials x 3 + x 2 + 1 and x 4 + x 1 + 1, corresponding to indices 5 and 6.

Fig. 3 .

 3 Fig. 3.34.: Sobol sets with polynomials of different degrees

Fig. 3 .

 3 Fig. 3.35.: Application of Gray Code on a Sobol sequence.

Fig. 3 .

 3 Fig. 3.36.: Hammersley point set [Ham60].

Fig. 3 .Fig. 3 .

 33 Fig. 3.37.: Overview of common state of the art samplers regarding spectral quality and discrepancy.

Fig. 4 . 4 . 1

 441 Fig. 4.1.: Owen's scrambling algorithm in 1-D [Owe95]

Fig. 4 . 3 .

 43 Fig. 4.3.: Owen's scrambling [Owe95] applied on a Sobol sequence [Sob67]. The measurements were done on a single realisation of this scrambling.

FrequencyFig. 4 . 4 .:

 44 Fig. 4.4.: Spectral comparison between Owen's scrambling, and our adaptation; both applied on the same Sobol Sequence. We can see that for small m values, artefacts still appear but as m increases, there is no longer any visible difference between our set and the original Owen.

Fig. 4 . 5 .:

 45 Fig. 4.5.: The measured discrepancy when changing the position of the red sample.

Figure 4 .

 4 Figure 4.5 shows how changing the position of a single sample affects the overall discrepancy.The discrepancy of the point set on the left is computed for various displacements of the red sample, and the results are plotted on the right. We can see that the lowest discrepancy is achieved as the red sample moves towards the origin, as intuitively expected when seeing the point set.

Figure 4 .

 4 Figure 4.6 on the next page shows the Fourier spectrum of a Blue Noise set obtained from the BNOT algorithm[START_REF] De Goes | Blue noise through optimal transport[END_REF], and the Fourier spectrum of this same set after optimization. We can note that a black cross appears in the middle of the spectrum. This cross appears as, to reduce the discrepancy, the samples aligned on either axis were pushed away from each other by our optimization.

Figure 4 .

 4 Figure 4.7 on page 86 presents the discrepancy of the original BNOT sampler, compared with its discrepancy after minimizing its generalized l 2 discrepancy. This Figure presents two graphs, showing the evolution of the l ∞ star discrepancy, and the evolution of the generalized l 2 star discrepancy.

Fig. 4 . 6 .Fig. 4 . 7 .:

 4647 Fig. 4.6.: Top: The original BNOT set, Bottom: Results of our discrepancy optimization. The black cross in the spectrum is characteristic in sets with lower discrepancy.

Fig. 4 .

 4 Fig. 4.8.: (0, k, 2)-nets with maximized minimum distance [Grü+08].

Algorithm 2 :Fig. 4 .Fig. 4 .

 244 Fig. 4.9.: Comparison of the results our genetic algorithm with a point set of 16 samples and an initial population of 200 random individuals, with the results from the exhaustive search.

Fig. 4 .

 4 Fig. 4.11.: Characteristics of our genetic algorithm.

Fig. 5 .

 5 Fig. 5.1.: Illustration of Equation 5.1 on the previous page, representing a stratified point set.

. 4 )

 4 A and B be 1-D disjoint subsets of [0, 1), thend(P n , A ∪ B) ≤ d(P n , A) + d(P n , B) .(5Note that the discrepancy measure d is here used for 1-D sets but is also defined in n-D.a := {a i } is a one-dimensional sequence of points in [0, 1). a(k) denotes the set of the first k elements of the sequence (b := {b i } and b(k), respectively). We consider a collection of k such sequences a (X) := a (X) i

Fig. 5 . 2 .:

 52 Fig. 5.2.: Notations for Lemma 1 on the preceding page.

5. 1

 1 Initial Stratified Low Discrepancy set allows us to call our construction a low discrepancy set. Note that our construction is only valid for 2-D point sets. We can thus construct 2-D low discrepancy stratified point set by picking any 1-D low discrepancy sequence S, and use it to supply the horizontal offsets along the columns: a X,Y := S Y , and the vertical offsets along the rows: b X,Y := S X .

Fig. 5 . 3 .:Fig. 5 . 4 .:

 5354 Fig. 5.3.: Our template 2-D stratified low discrepancy point set.

  .18)where each σ j is a permutation over {0...m -1}. (Fig.5.4) We know that this permutation have little impact on the discrepancy. Any subset of S differs from the corresponding subset of S in at most m points. According to [KN74, Theorem 4.1] (presented Section 4.1.2 on page 75), the difference in discrepancy of any subsets of length n > m of the two sequences is then bounded by m/n. Since m is fixed, the difference|D * (S ) -D * (S)|in discrepancy between the two sequences is O(1/n). Consequently, if S is an LD sequence, with O(log(n)/n) discrepancy bound, then the shuffled sequence S is also an LD sequence. The parameter m only alters the constant in the O(•) notation, but we want to keep m small so that we stay close to the discrepancy of S.

Fig. 5 . 5 .

 55 Fig. 5.5.: A 2-D stratified low discrepancy point set obtained from random scrambling of P φ .

5 . 5 )Algorithm 3 : 2 for each chunk j do 3 Y =j•m ; 4 Write the permuted vector back into the set; 5 end 6 end 7

 553234567 closely resembles a scrambled LD set like Owen's[START_REF] Owen | Randomly permuted (t, m, s)-nets and (t, s)-sequences[END_REF] (Fig.4.3 on page 79). Our localized permutation, however, admits much more control over the placement of individual points than Owen's axis-wise scrambling, and it enables real optimization. Algorithm 3 sums up our random scrambling. Stochastic permutation of P S .1 for each column X do Randomly permute the vector {a X,Y } j•m+m-1 Repeat with the rows by exchanging the roles of a, X and b, Y , respectively;

Fig. 5 . 6 .:

 56 Fig. 5.6.: Optimizing the samples placement trough dimension-wise rearrangement, using m = 8

Fig. 5 . 7 .:Fig. 5 .

 575 Fig. 5.7.: Bitwise representation of Eq. (5.19) and Eq. (5.20)

Fig. 5 .Fig. 5 .

 55 Fig. 5.9.: Characteristics of the Low Discrepancy Blue Noise sampler with a BNOT [DG+12] target profile. t = 128, m = 16 [Ahm+16b].

7 10

 10 Write the permuted vector Q x back to the output set; Repeat with the rows by exchanging the roles of a, X and b, Y , respectively;

Fig. 6 .

 6 Fig. 6.1.: Illustration of our discrepancy preserving 2-D scrambling for spectrum optimization.

Fig. 6 . 2 .:

 62 Fig. 6.2.: Characteristics of our sequence using admissible permutations

Fig. 6 . 3 .:Fig. 6 . 4 .:

 6364 Fig. 6.3.: Characteristics of our sequence using ad-hoc permutations and K=4. Note how a hierarchical construction increases the width of the black cross in the spectrum, this is due to the fact that the unaligned samples from the previous levels are more and more spaced out as we refine.

Fig. 6 .Fig. 6 . 7 .:

 667 Fig. 6.6.: 4-D point set generated with our system, together with two 2-D projections and their Fourier power spectra. Note that subsets of 16 points that occupy the same strata in one projection (red subsets in {(x, y)} projections and green subsets in {(u, v)} projections), occupy different strata in the other projection.

Fig. 6 .Fig. 6 .Fig. 6 .

 666 Fig. 6.8.: Subdivision of the domain [0, 1) s with increasing λ.

Fig. 6 .

 6 Fig.6.12.: s = 2, K = 2, λ = 1. Indices of the samples in S λ and of the samples in {s} λ r .

Fig. 6 .

 6 Fig.6.13.: s = 2, K = 2, λ = 2. Each vector v(i) is made by xoring the samples s (i) ∈ S λ-1 and p (i) ∈ P λ-1

Fig. 6 .

 6 Fig. 6.14.: s = 2, K = 2, λ = 2. The set Q λ is made by xoring each samples of each set {s} λ r with the vector v λ r associated to the tile T λ r

Fig. 6 .

 6 Fig. 6.18.: Binary representation of (x, y)| a

Fig. 6 .

 6 Fig. 6.20.: Geometric illustration of Property 6

Fig. 7 .Fig. 7 . 2 .

 772 Fig. 7.1.: l ∞ Star Discrepancy of various 2-D Low Discrepancy samplers.

Fig. 7 .

 7 Fig. 7.3.: l ∞ Star Discrepancy of various 2-D samplers.

Fig. 7 . 4 .:Fig. 7 . 5 .:

 7475 Fig. 7.4.: Comparing the l ∞ Star Discrepancy of our samplers with reference 2-D samplers.

Fig. 7 . 7 .: l 2

 772 Fig. 7.7.: l 2 Star Discrepancy of various 2-D Blue Noise samplers.

Fig. 7 .

 7 Fig. 7.8.: l 2 Star Discrepancy of various 2-D samplers.

Fig. 7 .Fig. 7 .

 77 Fig. 7.9.: Comparing the l 2 Star Discrepancy of our samplers with reference 2-D samplers.

Fig. 7 .

 7 Fig. 7.11.: Integration MSE and Variance measured on a Disk with the Sobol and Stratified samplers, with and without Cranley-Patterson shift.

Fig. 7 .

 7 Fig. 7.12.: Cornered Gaussian function without and with shifting. Discontinuities appear on the shifted function that did not exist in the original function.

Fig. 7 .

 7 Fig. 7.13.: Integration MSE and Variance measured on a Cornered Gaussian with the Sobol sampler, with and without Cranley-Patterson shift.

Fig. 7 .

 7 Fig. 7.14.: Variance of Monte Carlo estimators from Blue Noise samplers when integrating an analytical 2-D disk.

Fig. 7 .

 7 Fig. 7.15.: Variance of Monte Carlo estimators from Low Discrepancy samplers when integrating an analytical 2-D disk.

Fig. 7 .Fig. 7 .

 77 Fig. 7.16.: Variance of Monte Carlo estimators from various samplers when integrating an analytical 2-D disk.

Fig. 7 .

 7 Fig. 7.18.: Comparing the variance of Monte Carlo estimators of our samplers with reference samplers when integrating an analytical 4-D sphere.

  Results from a HDR Image showing a natural scene, with an O(ρ -1.59 ) profile.

Fig. 7 .

 7 Fig. 7.19.: Spectral characteristics of HDR integrands.

Fig. 7 .

 7 Fig. 7.20.: Variance of Monte Carlo estimators from Low Discrepancy samplers when integrating the inset hdr image.

Fig. 7 .Fig. 7 .

 77 Fig. 7.21.: Variance of Monte Carlo estimators from Blue Noise samplers when integrating the inset hdr image.

Fig. 7 .Fig. 7 .

 77 Fig. 7.23.: Comparing the variance of Monte Carlo estimators of our samplers with reference samplers when integrating the inset hdr image.

Fig. 7 .Fig. 7 .

 77 Fig. 7.25.: Variance of Monte Carlo estimators from Blue Noise samplers when integrating the inset hdr image.

Fig. 7 .

 7 Fig. 7.27.: Comparing the variance of Monte Carlo estimators of our samplers with reference samplers when integrating the inset hdr image.

Fig. 7 .

 7 Fig. 7.28.: Spectral properties of complex analytical arrangements. It is clearly visible in log scale that compared to the natural image from Figure 7.19, we still lack high frequencies.

  Properties of our random triangle arrangements O(ρ -2.25 ) profile.Properties of our aligned box arrangements O(ρ -2.34 ) profile.

Fig. 7 .

 7 Fig. 7.29.: Spectral properties of simpler analytical arrangements.

Fig. 7 .

 7 Fig. 7.30.: Comparing the MSE of Monte Carlo estimators of our samplers with reference samplers when integrating synthetic box arrangements.

Fig. 7 .

 7 Fig. 7.31.: Comparing the MSE of Monte Carlo estimators of our samplers with reference samplers when integrating synthetic triangle arrangements.

Fig. 7 .

 7 Fig. 7.32.: Ground truth reconstruction of the Zoneplate function.

Fig. 7 .

 7 Fig. 7.33.: Zoneplates of various Blue Noise samplers, using 262144 samples to reconstruct a 256x256 image.

Fig. 7 .Fig. 7 .Fig. 7 .

 777 Fig. 7.34.: Zoneplates of various Low Discrepancy samplers, using 262144 samples to reconstruct a 256x256 image.

Fig. 7 .Fig. 7 .Fig. 7 .

 777 Fig. 7.37.: Direct lighting (single light bounce) with 16 samples per pixel. The deterministic Sobolsampler is used on every dimension apart from the light plane. This dimension was sampled using various samplers, allowing to compare the impact of each of them and compare it with our sampler.

III 1 TFig. A. 5 .:

 15 Fig. A.5.: Discretizing a signal in frequency space

Fig. A. 6 .:

 6 Fig. A.6.: Reconstructing a discretized signal in frequency space

Fig. A. 7 .:

 7 Fig. A.7.: Reconstructing a discretized signal in spatial space

(A. 9 )

 9 Fig. A.9.: Usual filters in their frequency representation

  Fig. B.1.: Integral invariant computation (left) and notations (right) in dimension 2 [Coe+14].

  ball radius h in R = kh 1 3

  Fig. B.2.: Morton codes associated to cells of a linear quatree. Each morton code of a child cell is obtained by adding a suffix to its parent code (left). The adaptive representation consists of quadtree cells whose depth is view point dependent (middle). Finally, adaptive Marching Cubes is used to generate the triangulation (right).

  [START_REF] Gargantini | An effective way to represent quadtrees[END_REF]. This structure indexes each cell by a Morton code: the code of children cells are defined by the code of the parent suffixed by two bits (in dimension 2, three bits in dimension 3) (see Figure B.2-left).

Figure B. 2

 2 on the preceding page-right). Many authors have addressed this problem both for primal and dual meshes [Shu+95; SW04; LO10; Lew+10; Lob+14].

  Figure B.3 on the facing page illustrates our level of details (LoD for short) criterion. In dimension 2, if α denotes the viewing angle, an object at a distance d from the camera has a projected size on screen of 2 • d • tan(α). (see Figure B.3 on the next page-left).

Fig. B. 3 .

 3 Fig. B.3.: Notations (left) and adaptive meshing in dimension 2 using the LoD distance and angular criterion (right).

Fig. B. 4 .

 4 Fig. B.4.: GPU pipeline summary. Data buffers are represented in green and computations in red.Each computation retrieves data from a buffer and fills a new one.

  3 on page 175 is consistent with the triangulated surface obtained on GPU: triangles are used for visualization purposes but all computations are performed on the digital object G h (X) for estimators ( B.1 on page 176), ( B.2 on page 176), ( B.3 on page 176) and ( B.4 on page 176).To implement the integration on B R (x) ∩ X (Fig.B.1 on page 176), several strategies have been evaluated.

Fig. B. 5 .:

 5 Fig. B.5.: Here are the three approaches we used to evaluate the integrand. (left) is the naive approach,where we sum all the digital points that are inside the integrand domain (middle) shows the hierarchical decomposition of the ball, thus limiting the number of texture probes to do.(right) illustrates the dynamic approach, we probe the textures at a higher resolution and we refine it at each frame until we reach the finer level.

B. 6

 6 Figure B.6 on the next page-top shows re-Experiments

Fig. B. 6 .:

 6 Fig. B.6.: First row: Mean and Gaussian curvature estimation, (zoom of) first and second principal directions estimation on "OctaFlower" with a digital domain of 130 3 . Second row: Mean curvature computed in real-time (around 20 FPS) on a dynamic object.

  Figure B.6-bottom shows mean curvature estimation in real-time on a dynamic object. For this case, we simply evaluate the implicit expression at each vertex of a cell to decide if such cell is included or not into B R (x) ∩ X instead of updating the mipmap texture of densities at each time step. For all illustrations, we use directly normal vectors computed by algorithm discussed in Section B.3 on page 175.

H

  Fig. B.7.: Illustration of mean curvature computation on "OctaFlower" (digital domain of 130 3 ) using mipmap approximation with different levels: l 3 , l 2 , l 1 and l 0 (i.e. no approximation).

Fig. B. 8 .:

 8 Fig. B.8.: Timings in milliseconds (in logscale) obtained while visualizing an adaptive triangulation on three objects -"OctaFlower" (digital domain of 130 3 ), "Snow microstructures" (233 3 ) and "XYZ-Dragon" (510 3 ) -by computing the curvature with a regular grid (left) and with hierarchical algorithm (right), with two different radii for each object. In orange color: time required to extract the triangulation. In blue color: time required to compute the curvature tensor at different levels: from l 4 (light blue) to l 0 (dark blue). Timings are given using a NVIDIA GeForce GTX 850M GPU.Tutu

  3 on page 35), which will let us underline where do our contributions fit in (Chapter 4 on page 74, Chapter 5 on page 93, and Chapter 6 on page 107). Finally, we will present how those samplers behave in practice through experimental results (Chapter 7 on page 133). Note that is this PhD, we study both sequential and non-sequential samplers. However, only infinite and stationary point processes are considered, with most of them being also isotropic.

	http://cpc.cx/jDK
	http://cpc.cx/jDL
	http://cpc.cx/jDM
	I'm David Coeurjolly and I approve this manuscript
	Chapter 1
	s

  . It uses either Lloyd's relaxation or a Dart Throwing algorithm to generate a sampling pattern matching a particular Pair correlation function (Section 2.3 on page 25). For the Dart Throwing variant, the method changes the conditions of acceptance or rejection of the samples. When adding a new sample, the PCF (see Section 2.3 on page 25) of the new set is computed. The distance between this PCF and the targeted PCF is computed using a l ∞ or a l p distance between the two PCFs. A sample is accepted if this distance between PCFs is below a certain threshold. Similarly,[START_REF] Öztireli | Analysis and Synthesis of Point Distributions based on Pair Correlation[END_REF] adapted the Lloyd's algorithm by weighting the gradient in the optimization with the PCF. As it is based on simple Dart Throwing or simple gradient descent, this method cannot generate massive amounts of samples. Furthermore, this method finds it more difficult to generate very regular structures, as the PCF (when used with Euclidean distance) is an isotropic measure. But, as the PCF is also a s-D measure, this methods can generate pointsets in higher dimensions.

  The a k values are defined in the F 2 domain. This domain is a ring over Z/2Z with operators for addition ⊕ and multiplication ⊗. The C d matrices, defined in F N ×N 2 , are called generating matrices. From here, what differentiates most methods generating digital nets and digital sequences are the C d matrices used.

  Chapter 5 on page 93 will now present a first Low Discrepancy Blue Noise sampler (LDBN) that was developed jointly by us and the Universities of Konstanz and Shengen. This sampler is a first step towards hybrid samplers but it is still limited to 2-D and does not present any sequential properties. Then, Chapter 4 on page 74 will present simplistic methods (inspired from energy minimization or genetic algorithms) to develop multi-dimensional low discrepancy sequences with controllable spectrum. Finally, Chapter 6 on page 107 will improve upon those naive ideas to present a new scrambling technique applied on the Sobol sequence that preserves its sequentiality while both giving a Blue Noise spectrum to its 2-D projections and maintaining their Low Discrepancy property. It also maintains some uniformity of the s-D set.

		tutu
	Good
	Spectrum	BNOT
		[DG+12]
		Canonical
		Blue Noise
	Stratified	Poisson Disk
		Low Discrepancy
	Whitenoise	Dart Throwing
		Halton
		[Hal64]
		Sobol
	Bad	[Sob67]
	Spectrum
	High	Low
	discrepancy	discrepancy

  Size of our lookup table in number of entries need to generate a certain amount of samples, depending on the K factor used. Please note that this is in logarithmic scale.

		10000
								K=2
								K=4
								K=8
		1000
	Number of patterns needed	100
		10
		1
		1	10	100	1000	10000 100000 1e+06	1e+07	1e+08	1e+09
					Number of samples generated
	Fig. 6.5.:

  We start byproperly defining what is a continuous signal, and how it differs from a discrete signal. We denote continuous signals as continuous-time (CT) signals. Those signals (despite their name, which stands mostly for historical reasons) do not necessarily have anything to do with time. They are the family of signals defined from analytic functions; if f is a continuous signal, then an infinitesimal difference in x leads to an infinitesimal difference in f (x). Note that those signals are not necessarily differentiable everywhere. Those signals are often used to modelize the real world and correspond to theoretical definitions. As their domain of definition is infinite, as soon as you try to evaluate such signals with a computer, you have to turn then into discrete-time (DT) signals.A DT signal can be seen a set of s-dimensional samples with their associated values. It is not continuous, not defined on a infinite domain, and not differentiable. Since a computer cannot represent the infinite amount of existing numbers (due to the floating precision limit), they are the only kind of signals manageable on a computer. Numerizing a CT signal implies that this signal will be discretized over a set of computer-representable values. Even though most often this is ignored as this floating point precision is higher than what we consider negligible error, we will see that when handling signals with an infinite maximal frequency, this becomes a problem. The difference between CT and DT signals is illustrated Figure A.1.

	f	III T	f (t) • III T	t
	Fig. A.4.: Discretizing a signal in spatial space
	ω
	(a)	t	(b)	t

Fig. A.1.: (a) A 1-D CT signal f (t) and (b) its equivalent DT signal ḟ (t) t × t =

  We call this limit the Nyquist limit. If a function has been discretized with T < 2ω, it is said to be undersampled. It is easy to see that in frequency space, this causes the interval of the frequency Dirac comb to be too small; the aliases will cover each other, creating what are called aliasing artefacts. Aliasing might also occur on properly sampled signals, if the reconstruction process fails to capture correctly the Fourier spectrum of the original function (Fig. A.10).All this theory allows proper reconstruction of a 2-D image from a set of discrete radiance values. However, computing those radiance values in itself is also a challenge, as detailed in Chapter 1 on page 1.

	Interactive Curvature Tensor	B
	Visualization
	B.1 P.h.D. Context
	T ≥ 2ω.	(A.10)

  In this topic, many solutions have been developed in Computer Graphics [Shu+95; SW04; LO10; Lew+10; Lob+14]. The method developed by Lengyel et al. [LO10] suits best our needs. It combines an octree space partitioning with new Marching Cubes configurations to generate adaptive meshes on CPU from static or near static data. We propose here a full GPU pipeline inspired by this method, that maintains a view dependent triangulation. Our high framerate allows us to inspect dynamic 3D data in real-time.

  Such curvature estimators have multigrid convergence properties[START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF]: when the digital object becomes finer and finer, i.e. when the digitization step h tends to zero, the estimated quantities on ∂[G h (X)] h converges (theoretically and experimentally) to the associated one on ∂X in O h

	κR (Z, x, h) :=	3π 2R	-	3 Area(B R/h (x/h) ∩ Z, h) R 3	,	(B.1)
	ĤR (Z, x, h) :=	8 3R	-	4 Vol(B R/h (x/h) ∩ Z, h) πR 4	.		(B.2)
	1
	3

where B R/h (x/h) is the ball of digital radius R/h centered on digital point (x/h) ∈ Z.

  Comparison of number of texture fetches, L 2 and L ∞ error obtained on "OctaFlower" with a digital domain of 130 3 when computing mean curvature with two radii: 8 and 16, with hierarchical algorithm (H) and l ≥ l 0 approximation algorithms. The object is triangulated at full resolution with a regular grid and contains 282,396 vertices.

				l 4	l 3	l 2	l 1	l 0
	Number of	R = 8	1468		2	28	260 2104
	texture	R = 16	5706	2	28	90	2120 17080
	fetches L ∞ error	R = 8	0.051		0.306	0.085	0.047	0
	(w.r.t. l 0 )	R = 16	0.053	0.146	0.041	0.017	0.005	0
	L 2 error	R = 8 3.51e-05		2.67e-04 7.57e-05 2.02e-05	0
	(w.r.t. l 0 )	R = 16 4.43e-05 1.39e-04 4.16e-05 1.57e-05 2.07e-06	0
	Tab. B.1.:

  Présente l'ensemble de points 2-D initial de notre méthode. Cet ensemble est permuté sur chaque coordonnées (b,c) Pour donner un ensemble qui soit bruit bleu à l'arrivée (d). On peut noter que l'ensemble de départ a un spectre de Fourier très différent de celui de l'ensemble d'arrivée. On peut également noter que le spectre de Fourier de l'ensemble d'arrivée présente une croix en son centre. Cette croix est symptomatique de la basse discrépance et n'est en aucun cas un handicap aux capacités de notre échantillonneur en terme d'intégration. . . . . . . . . . . . . . . vii 0.2 (a) Présente l'ensemble de points 2-D initial, répartis dans des tuiles à un niveau l. On permute individuellement chaque tuile (b,c) pour obtenir un ensemble qui soit bruit bleu à l'arrivée (d). Une tuile contient toujours le même nombre de points. Pour rajouter des points, on subdivise les tuiles. Cette subdivision peut être faite localement pour faire varier la densité de points (e). . . . . . . . . . . . . . . . . . viii 1.1 The rendering process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Different artefacts observed when rendering a scene using different sampling strategies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Various 2-D point sets P 1024 generated from different samplers S. . . . . . . . . . 2.1 Spectral test in 2-D on Korobov lattices. Figure from [Lem09] . . . . . . . . . . . 2.2 Congruent boxes in 2-D on the Sobol sequence [Sob67] . . . . . . . . . . . . . . . 2.3 The CDF of a pointset P n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The Uniform distribution U(a, b) . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 The Kolmogorov-Smirnov difference (in red) between the CDF of a point set (in blue) and the CDF of the uniform distribution (in black). . . . . . . . . . . . . . . 2.6 Different measures for discrepancy in 2-D. . . . . . . . . . . . . . . . . . . . . . . 2.7 Computing the discrete l ∞ star discrepancy in 2-D. . . . . . . . . . . . . . . . . . 2.8 Fourier Serie decompositions F N of the square wave function and their corresponding Fourier transforms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Power spectrum and Radial average characterization. . . . . . . . . . . . . . . . . 2.10 Illustration of Ripley's K function and of the PCF . . . . . . . . . . . . . . . . . . 2.11 Ripley's K-function, the L-function and the PCF of the pointset P n . . . . . . . . . 2.12 Product of the radial spectra of the sampling pattern Ps and of the integrand PI . 2.13 Blue Noise fourier spectrum and radial average. . . . . . . . . . . . . . . . . . . . 2.14 Comparison of convergence speeds. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Characteristics of a Whitenoise sampler. . . . . . . . . . . . . . . . . . . . . . . . 3.2 Comparison between a Whitenoise point set and a Stratified point set. . . . . . . 3.3 Characteristics of a Stratified sampler. . . . . . . . . . . . . . . . . . . . . . . . .

	0.1	(a)

Nous avons pu produire deux travaux sur ce domaine. Le premier présente un premier échantillonneur en 2-D qui ne présente pas d'aliassage tout en ayant une faible erreur lors de l'approximation de l'intégration. Il est décrit Section 0.3. Section 0.4 présente notre deuxième contribution. Nous y présentons un échantillonneur toujours à faible variance et anti-aliassé, mais pour de plus hautes dimensions. Il permet également de contrôler localement la quantité d'échantillons et offre la possibilité d'augmenter incrementalement le nombre d'échantillons.

1.1 The Rendering process

1.3 Contributions

1.4 Basic definitions

2.1 Uniformity

2.2 Spectral analysis

-A purely random pointset (right), where every frequency appears. Its power spectrum (middle) is isotropic as all frequency appear equally, and its radial average (left) confirms the absence of peaks.6 -

2.3 Spatial Analysis

2.4 Relationship to error in integration

2.5 Conclusion

Chapter 2 Characterizing a set of samples

3.1 Historical samplers

3.2 Poisson Disk and Blue Noise samplers

3.3 Low discrepancy samplers

4.3 From an Anti-aliased to a Low Discrepancy set

Auto-differentiation tools use operators overloading to obtain the derivative h (y) at a point y of a code expression h := f (g(k(x))).Basically, the tool uses an automatic chain rule to compute the derivative at y back-propagating the values through the source code graph (https://en.wikipedia.org/wiki/Automatic_differentiation) 4.3 From an Anti-aliased to a Low Discrepancy set

4.4 From a Low Discrepancy to an Anti-aliased set

4.5 Conclusion

5.2 Scrambling

6.1 Overview

6.2 Technical description

6.4 Conclusion

Chapter 6 Multi-Dimensional Sequences with Low-Discrepancy Blue-Noise 2-D Projections

7.1 Discrepancy

7.2 Integration error

7.3 Aliasing

Chapter 7 Experiments

A.1 Discretization

A.3 Reconstruction

B.3 Curvature Tensor Estimation

B.4 Isosurface Extraction on GPU

B.5 Interactive Curvature Computation on GPU

Chapter B Interactive Curvature Tensor Visualization

B.7 Conclusion and Discussion

Following this, Chapter 6 on page 107 present a new construction to generate sequentially low discrepancy sets in higher dimensions.

Conclusion

From Eq ( 6.3 on page 120), we have

where v λ r is a singleton defined as v λ r = {s} λ-1,λ r ⊕ {p} λ-1,λ r (6.21) with {s} λ-1,λ r the i th sample in S λ-1 and {p} λ-1,λ r the i th sample in P λ-1 .

Therefore

Then, from equation 6.4 on page 122 we have

and by property 4 on page 122, we know that Π preserves the position of the samples from P λ-1 . Therefore, is P λ-1 ⊂ Q λ , P λ-1 ⊂ P λ .

Our local permutations are admissible permutations

In this section, we demonstrate that the set of permutation {π λ r } built by our method (See Section 6.1.2 on page 110) is a subset of the admissible set of permutation Π.

Lemma 7

The local permutations {π λ r } as defined in Section 4.1 of the main article define an admissible Π.

Proof 6

We already explained how our permutation can preserve the position of a given marked sample (constraint (i) leading to Lemma 6 on the preceding page). We give here more details on how we exchange the trailing bits, to demonstrate that we preserve the net properties of the initial set and its 1-D projections. A DT signal is obtained from a CT signal by computing the product of the CT signal with particular signals, namely the Dirac impulse and the Dirac comb.

Tutu t

A.1 Discretization

The Dirac impulse (Fig. A.2) is a CT signal defined as

The issue here is that δ(0) has therefore an undefined value. Since the point 0 is infinitely small, to ensure both Equation ( A.1) and Equation ( A.2), we need to have δ(0) such that δ(0) ∞ = 1 which is undefined. To solve this issue, the Dirac impulse can be redefined as a DT signal with A discrete signal can also be reconstructed back to a CT signal, but only under some conditions.

To understand those conditions, we need to look at what happens in the frequency space when we discretize a signal.

6.6

4-D point set generated with our system, together with two 2-D projections and their Fourier power spectra. Note that subsets of 16 points that occupy the same strata in one projection (red subsets in {(x, y)} projections and green subsets in {(u, v)} projections), occupy different strata in the other projection. . . . . . . . .

6.7

Comparison between an original 8-D Sobol set and this point set scrambled using our method. Note that we used K = 4, and thus can only apply our scrambling to the the indices (1, 2) and (3,7). The remaining dimensions were permuted using Owen', as described in Chapter 4 on page 74. . . . . . . . . . . . . . . . . . . . .

6.8

Subdivision of the domain [0, 1) s with increasing λ. . . . . . . . . . . . . . . . . .

B.5

Here are the three approaches we used to evaluate the integrand. (left) is the naive approach, where we sum all the digital points that are inside the integrand domain (middle) shows the hierarchical decomposition of the ball, thus limiting the number of texture probes to do. (right) illustrates the dynamic approach, we probe the textures at a higher resolution and we refine it at each frame until we reach the finer level.

List of Tables

1.1

Table of the notations we will use throughout this document. . . . . . . . . . . . .

3.1

Worst case convergence speeds of the integration error when using warious sampling patterns, along with their discrepancy. . . .