École Doctorale

D E Mathématiques

POUR OBTENIR Arthur Dufay

Xavier Granier Professeur

. . . . . . . Professeur, DIRO Pierre Poulin

XLIM Daniel Meneveaux Professeur

Mathias Paulin

Professeur

Je Tiens Tout D'abord À Remercier

Daniel Méneveaux

Gaël

THÈSE PRÉSENTÉE À L'UNIVERSITÉ DE BORDEAUX

Keywords: Rendu, Carte Graphique, Rasterisation, Eclairage Global, Lancer de Rayons, Path Tracing Rendering, GPU, Rasterization, Global Illumination, Ray Tracing, Path Tracing

Titre Rendu adaptatif haute-qualité d'environnements virtuels à photométrie complexe Résumé La génération d'images de synthèse pour la production cinématographique n'a cessé d'évoluer durant ces dernières décennies. Pour le non-expert, il semble que les eets spéciaux aient atteint un niveau de réalisme ne pouvant être dépassé. Cependant, les logiciels mis à la disposition des artistes ont encore du progrès à accomplir. En eet, encore trop de temps est passé à attendre le résultat de longs calculs, notamment lors de la prévisualisation d'eets spéciaux. La lenteur ou la mauvaise qualité des logiciels de prévisualisation pose un réel problème aux artistes. Cependant, l'évolution des cartes graphiques ces dernières années laisse espérer une potentielle amélioration des performances de ces outils, notamment par la mise en place d'algorithmes hybrides rasterisation/lancer de rayons, tirant prot de la puissance de calcul de ces processeurs, et ce, grâce à leur architecture massivement parallèle.

. Pour nir, nous avons tenté de combiner l'algorithme de Path Tracing et les solutions Many Lights, toujours dans le but d'améliorer la prévisualisation de l'éclairage global. Cette thèse a aussi donné lieu à la soumission de trois mémoires d'invention et a permis le développement de deux outils logiciels présentés en Annexe A.

Introduction:

L'informatique graphique est la science qui regroupe tous les méthodes de communication visuelles via un ordinateur. Un de ces sous domaine, la synthèse d'image, ou rendu, et l'ensemble des techniques qui permet la création d'images de synthèse à l'aide d'un ordinateur. La synthèse d'image à plusieurs champs d'applications: les jeux vidéos, l'architecture, la publicité, les films d'animations, les effets spéciaux ... Au cours des dernières décennies, l'utilisation d'effets spéciaux dans les films, ainsi que la production de films d'animation n'a cessé de croître. Les effets spéciaux ont fait un énorme progrès depuis l'apparition du premier film utilisant de l'informatique graphique 3D:

Tron de Disney (1982) . Pour le non-expert, il semble que les effets spéciaux aient atteint un niveau de réalisme ne pouvant être dépassé. De nos jours, les images de synthèse se fondent naturellement dans le film. Les films d'animation ont aussi fait un grand bond en avant depuis Toy Story de Pixar (1995) . Bien que la puissance de calcul disponible pour générer ce type d'images ait drastiquement augmenté, la quantité de données requises pour générer des images avec un tel niveau de réalisme croît aussi rapidement. En plus de cela, les algorithmes et techniques impliqués dans la génération de ces films sont en constante évolution. C'est particulièrement vrai pour la prévisualisation d'effets spéciaux. Les artistes, et plus spécifiquement les lighters, sont encore limités par les outils logiciels qui leur sont fournis. A l'heure actuelle, les solutions logiciels disponibles pour la prévisualisation d'effets spéciaux sont soit trop lent, soit ne correspondent pas en terme d'apparence au images finales du film. Elles ne donnent pas un retour appréciable aux artistes.

Cette thèse tente d'améliorer les capacités d'un outil logiciel de production d'effets spéciaux développé à Technicolor. Pour cela, nous améliorons un moteur de rendu développé en interne: 3DCast , livrés à des entreprises productrices d'effets spéciaux telles que MPC (The Moving Picture Company) .

Organisation du manuscrit de thèse:

Ce manuscrit est divisé en cinq chapitres, chacun commençant par un rapide résumé de son contenu. Ils décrivent les différentes parties de la solution logicielle pour la prévisualisation d'effets spéciaux développées au cours de cette thèse. En effet, développer un moteur de rendu permettant le calcul d'illumination globale requiert de couvrir différents domaines d'études tels que: la représentation des matériaux, le rendu, l'intégration de Monte-Carlo, les structures de données d'accélération spatiale et bien d'autres.

Résumé des contributions:

Premièrement un moteur de rendu se basant sur l'algorithme de tracé de chemins à été implémenté au sein de la plateforme 3DCast . Avec, une solution pour améliorer l'interactivité se basant sur un quad-tree en espace image (cf. Section 2.2.4).

Un brevet sur l'amélioration du tracé de rayons d'ombrage a été soumis (cf. Section 3.6).Il groupe les rayons d'ombrage par sources de lumières et lance les rayons de façon inverse (des sources de lumière vers les surfaces) pour améliorer la performance du tracé de rayon.

Une optimisation de l'algorithme de traversée de la hiérarchie de boîtes englobantes sur carte graphique a été proposé dans la Section 4.3.2. Elle utilise un encodage de la traversée au sein de la structure accélératrice et permet d'accroître l'efficacité de l'algorithme de tracé de rayons dans notre moteur de rendu.

Cette thèse a aussi permis le développement d'une nouvelle méthode de décorrélation des échantillons qui s'est prouvée être avantageuse pour plusieurs domaines d'intégration multidimensionnelle incluant le tracé de chemins et le calcul de SSAO. Cette contribution a été publiée dans une conférence majeure et soumis comme brevet.

Finalement, pour améliorer le calcul d'ombres, un brevet permettant d'hybrider shadow-maps et lancer de rayons a été soumis, il est présenté en Annexe B de ce manuscrit.

Travaux futurs:

Travaux de développement:

Plusieurs sous partie de notre moteur de rendu non pas pu être achevées pendant cette thèse. Cela mène à divers travaux futurs de développement.

Premièrement le support de matériaux plus complexes, par exemple les BSSRDF (rendu de peau), est essentiel pour le rendu de surface complexe rencontrable dans un cas de production cinématographique. De pair avec ces matériaux, nous voudrions aussi implémenter un autre type de pipeline de rendu, en suivant la philosophie de wavefront proposé par [START_REF] Laine | Megakernels considered harmful: Wavefront path tracing on gpus[END_REF]. Cela serait particulièrement profitable avec des matériaux plus complexes.

Pour la structure d'accélération spatiale nous aimerions nous intéresser aux structures plus complexe, notamment en ajout la possibilité de découper les triangles à notre hiérarchie de volumes englobants. Cela permettrait d'améliorer les performances de celle ci et lui permettrait de se rapprocher des performances obtenus avec un kd-tree.

Finalement, bien que la flexibilité et la rapide mise en place proposé par les Compute Shaders d'OpenGL nous ai convenu, nous aimerions implémenter une version se basant sur le framework NVIDIA CUDA. Nous pensons que cela permettrait de plus finement pouvoir déboguer et profiler notre code grâce à tous les outils fournis par ce framework.

Travaux de recherche:

L'axe majeur de recherche sur lequel nous aimerions nous concentrer est les solutions bidirectionnelles tels le tracé de chemin bidirectionnel, ou un une solution plus complexe, par exemple en combinant tracé de chemin et tracé de chemin bidirectionnel. Nous pensons qu'en utilisant des heuristiques pour décider si un chemin caméra doit se connecter à un chemin lumineux nous pouvons améliorer la solution. Une de ces heuristiques pourrait être par exemple de se baser sur les matériaux rencontrés le long du chemin. Un chemin spéculaire peut être difficile voir impossible à connecter avec une autre BRDF spéculaire. Une autre heuristique pourrait se baser sur la quantité d'énergie lumineuse transportée par le rayon. Par exemple ne connecter que les chemins caméra ne transportant peu d'énergie à des sources de lumières. Utiliser un tel algorithme pourrait grandement améliorer la performance d'un moteur de rendu basé sur l'algorithme de tracé de chemin bidirectionnel en permettant d'avoir des chemins unidirectionnels et d'autres bidirectionnels. En effet, en couplant les deux nous pourrions gagner en temps de calcul.

Introduction Motivations

Computer Graphics (CG) is the science that regroups all the visual communication methods via a computer. One of its elds, image synthesis, or rendering, is the set of techniques that enable the creation of synthetic images with a computer. Image synthesis has many applications: video games, architecture, advertisement, computer-animated movies, Visual Eects (VFX) ... (cf. Figure 1). Over the last decades, the use of Visual Eects in movies, as well as the production of computer-animated movies have become more and more prominent. VFX have made tremendous progress since one of the rst movies featuring 3D CG: Disney 's Tron (1982). To the eye of a non-expert, we seem to have reached a level of realism that cannot be outperformed. CG images now blend seamlessly into lmed images (cf. Figure 2). 3D computer-animated lm also made impressive progress since one of Pixar's Toy Story (1995) (cf. Figure 3). Even though the computing power available to generate such movies has drastically increased, the required amount of data to render images with such realism is also growing fast. In addition to that, the algorithms and techniques involved in the making of such movies are still evolving. This is especially true for VFX previsualization. VFX artists, especially lighters, are still limited in their work by the tools we provide them. At this time, the available software Figure 2: Left: CG at its early stages in one of the rst movies featuring 3D CG: Disney 's Tron (1982). Right: Impressive realistic CG integration in the recent movie The Martian (2015). Top: green chroma-key without VFX, Bottom: insertion of VFX. solutions for VFX previsualization (previz) are either too slow or do not match accurately the look of the nal movie images. They do not give an appreciable feedback to artists. This thesis tries to enhance the capability of VFX production tools developed at Technicolor. To do so, we improve an in-house rendering engine: 3DCast (cf. Chapter 2), delivered to VFX companies such as The Moving Picture Company (MPC).

deed, developing a rendering engine with global illumination features requires covering various elds of study such as material repsesentation, rendering, Monte-Carlo integration, spatial acceleration data structures, and so on.

Chapter 1 is dedicated to the necessary theoretical background. It presents the dierent concepts of Computer Graphics used during this thesis such as rendering, rasterization, ray tracing, and global illumination techniques. A mathematical background on Monte-Carlo integration is also presented. This chapter explains, after reviewing dierent global illumination techniques, why we chose to focus our eort on path tracing during this thesis.

Chapter 2 presents the industrial context of this thesis. It also introduces Computer Graphics notions such as materials and light sources used in this thesis. In addition, it explains the architecture of the path tracer we developed, as well as some techniques we implemented to accelerate it.

Chapter 3 focuses on GPU architecture and its use. In this chapter, we get further into details of the GPU implementation of our path tracer. First, we describe in detail the GPU architecture in order to give a clear view of its potential and limitations. Then, two implementation schemes are presented. Finally, we describe a patent we submitted, to accelerate shadow ray queries in rendering pipelines using ray tracing.

Chapter 4 deals with another crucial point of a rendering engine: spatial acceleration data structures. Several of them are presented to motivate further our choices. It also presents the more advanced techniques we implemented using a BVH on the GPU.

Chapter 5 is dedicated to our micro jittering technique [START_REF] Dufay | Cache-friendly micro-jittered sampling[END_REF].

It starts by explaining the challenges of random number generation and their uses on GPU. After reviewing quickly the caveats of dierent random number sequences, we present our new GPU cache friendly decorrelation technique.

All the contributions, publications, and patents published are summarized in the Conclusion of this thesis. This manuscript ends by discussing some potential future work.

Chapter 1 Theoretical Background

The essential part of this thesis relies on computing global illumination (GI) by solving the rendering equation. In this sense, we present in this chapter the required background to understand our work. We start, in Sections 1.1 and 1.2, by giving the denitions of rendering and global illumination. Then, in Section 1.3, we present the mathematical tools needed to understand our work. Sections 1.4 and 1.5 are dedicated to two main solutions to achieve rendering: rasterization and ray tracing. Finally, we get further in detail in Section 1.6, by reviewing the existing global illumination solutions that have been introduced in the last few years.

1.1. Rendering

Rendering

In Computer Graphics, we dene the process of "rendering" (cf. Figure 1.1) as the generation of an image (an array of pixels) from the description of a scene. A scene is described by four main components:

• a camera model: a virtual camera composed of a 3D position, a viewing direction, a eld of view and a resolution.

• some geometric data (e.g., mesh).

• a set of materials that control the scattering (reection and/or refraction) properties of the geometry.

• a set of light sources, which are described by geometrical and emissivity properties.

These geometric data could be for instance a set of 3D points, the denition of a volume, or in the context of this thesis, a set of 3D polygonal meshes. Rendering, also known as image synthesis, has many applications: video games, architecture, CG movies, Visual Eects (VFX), data visualization, Virtual Reality, ... Each eld of application has its own specic demands on features and image quality depending on the computing power capability and the required framerate. In this thesis, we focus on using rendering for previsualization of VFX for movie production. This implies complex lighting eects, that we describe later in this document, computed on massive 3D scenes containing high resolution 3D meshes (several millions of triangles).

Rendering can be achieved by several techniques that can be divided into two types of methods. The ones based on rasterization (cf. Section 1.5) and the ones based on ray tracing (cf. Section 1.4). We describe more in detail the pros and cons of these methods in the following sections of this chapter. Later, we will see that the two of them can be combined to take advantage of both. rithms aim at producing highest quality images intended to be integrated into a movie using a compositing stage. Computing images for nal render is mostly done on render farms. The latter are clusters of computers containing highend CPUs and huge amount of memory to avoid as much as possible swapping (transfer from hard drive to CPU memory) that really slows down the process of rendering by starving CPU threads. To this date, VFX studios still rely on CPU render farms and GPUs render farms exist but are still at their early stages. They are not used in production mostly because of their high power consumption, which drastically increases electricity costs. GPUs are also much more limited on memory than CPUs1 , which makes their use in a production renderer more tedious, even though out-of-core rendering algorithms exist. Rendering just one image, such as the one shown in Figure 1.2, on a render farm for a movie could take from hours to days, hence the need of a previsualization algorithm: it greatly helps VFX artists to have a real-time preview of the scene they are working on.

1.2. Global Illumination

Global Illumination

In this thesis, we address the problem of rendering 3D scenes with complex lighting eects (cf. Global illumination is crucial to obtain images that look realistic. We refer to direct lighting as light contribution that comes directly from a light source (see Figure 1.3), bounces once on an object of the scene before hitting the camera sensor. Indirect lighting represents light contribution from rays that hit more than one object, before hitting the camera sensor.

Path Classication

Light paths may be classied using a regular expression, introduced by Heckbert [1990]. It describes the materials encountered at each vertex of a light path using the following notations:

• L: a light source

Global Illumination

For instance, a path noted LS + DE refers to a path that starts on a light source, hits 1 to n specular surfaces, then a diuse surface, and nally the camera sensor. Such path represents a caustic, see Figure 1.5.

Radiometric Units

We introduce here some radiometric quantities and their I.S. units, needed for the comprehension of this document.

A Steradian noted sr is the unit of a solid angle, it is dened as the ratio between the area subtended and the square of its distance from the origin (see Figure 1.6). The Irradiance noted E is the radiant ux received by a surface per unit area.

Its unit is in

W • m -2 .
The Radiosity noted B is the radiant ux leaving a surface per unit area. Its

unit is in W • m -2 .
The Radiance noted L is the radiant ux emitted, reected, transmitted or received by a surface, per unit solid angle per unit projected area. Its unit is in W

• sr -1 • m -2 .

BRDF

In Computer Graphics the reective behavior of materials is described by a 4D function called Bidirectional Reectance Distribution Function (BRDF), introduced by [START_REF] Nicodemus | Directional reectance and emissivity of an opaque surface[END_REF], noted f r or brdf , and that is dened as:

f r (ω i , ω o ) = dL r (ω o ) dE i (ω i ) = dL r (ω o ) L i (ω i ) cos θ i dω i (1.1)
where ω i refers to the incoming light direction and ω o the outgoing or reected direction. The BRDF holds the ratio of reected radiance along ω o to the irradiance incident on the surface from direction ω i . Its units is in inverse 1.2. Global Illumination steradian (sr -1 ).

A physically based BRDF must respect three constraints:

• BRDF cannot create energy, meaning that Ω f r (ω i , ω o ) cos θ i dω i 1.

• it must be positive: f r (ω i , ω o ) 0.

• as stated in [START_REF] Stokes | On the perfect blackness of the central spot in Newton's rings, and on the verication of Fresnel's formulae for the intensities of reected and refracted rays[END_REF] and [START_REF] Helmholtz | The Helmholtz reciprocity principle[END_REF], for non magnetic materials, it obeys Helmholtz reciprocity:

f r (ω i , ω o ) = f r (ω o , ω i ).

The Rendering Equation

To compute global illumination, Kajiya introduced [START_REF] Kajiya | The rendering equation[END_REF] the rendering equation:

L o (x, ω o ) = L e (x, ω o ) + Ω f r (x, ω i , ω o )L i (x, ω i )(ω i • n) dω i (1.2)
which states that the outgoing radiance towards ω o at x on surface S is equal to the radiance emitted from S at x towards ω o plus the sum of all incoming light that is reected by S. To compute the global illumination we have to integrate all the incoming radiance over a hemisphere centered at x, oriented towards n (the normal of S at x) (see Figure 1.6).

This equation can be rewritten as

L o (x, ω o ) = L e (x, ω o ) + T (L o (x, ω o )) L o (x, ω o ) = L e (x, ω o ) + T (L e (x, ω o )) + T 2 (L o (x, ω o )) L o (x, ω o ) = ∞ i=0 T i (L e (x, ω o )) (1.3)
This reformulation using T , the light transport operator introduced by [START_REF] Veach | Robust Monte-Carlo Methods for Light Transport Simulation[END_REF], and n, the dimension of the ray space, exposes the high recursivity of the rendering equation. Since the dimension of the ray space is potentially innite, deterministic methods are not really suited for solving the rendering equation. In fact their convergence rate is in O(n -c d ) where c depends on the integration scheme and d is the dimension of the space of integration. Furthermore solving the rendering equation is not trivial since there is, in general, no analytical solution.

Computing this integral also involves solving the visibility problem between surfaces. Some work by [START_REF] Durand | The visibility skeleton: A powerful and ecient multi-purpose global visibility tool[END_REF] has been done to precompute visibility but it solves the problem for static geometry only. Furthermore, the complexity of their algorithm is in O(N 5 ) in time and O(N 4 ) in memory, N being the number of triangles in the scene, which is not acceptable in our context. To be fair, it still allows some possible editing of the scene: materials and light sources can be modied.

Theoretical Background

In addition to this integral, we also need to integrate over the area of a pixel in order to simulate the camera lens behavior and its complex eects such as depth of eld (DOF) or Bokeh.

When dealing with animated scenes, one also needs to integrate over the exposure time to simulate motion blur.

Finally, for more advanced rendering, it is also crucial to integrate over the spectrum of light, but this is not the subject of this thesis. Spectral rendering is more intended to architecture and material design.

For all these reasons, in general, Monte-Carlo integration is used.

Monte-Carlo Integration

Here, we present the mathematical background on Monte-Carlo techniques needed to understand this thesis.

Denitions: estimator, pdf and cdf Monte-Carlo integration is a method to integrate a function when no analytical solution exists. It relies on the ability to evaluate the function at random positions of its denition domain. It integrates it by summing up a set of samples of this function. This sum is written as follows:

I = Ω f (x)dx ≈ Q N Q N = 1 N N i=1 f (x i ) pdf (x i ) (1.4)
Q N is an estimator of the integrand I, by the law of large numbers we have

lim N →∞ Q N = I, (1.5)
which shows that the estimator converges to the correct solution. The pdf (x i )

is the probability density function, a function whose integrand is equal to the probability of choosing sample x i in the domain of integration Ω.

The cumulative distribution function of a random variable X is noted cdf (x)

and is dened as:

cdf (x) = P r{X x} (1.6)
The pdf and cdf are related by:

β α pdf (x)dx = P r{α X β} = cdf (β) -cdf (α) (1.7) or pdf (x) = dP r(x) dx (1.8) 1.3. Monte-Carlo Integration
Each sample x i is randomly chosen. For that several sampling methods exist. A sampling method has an associated pdf . In practice, the inverse cdf , noted cdf -1 , is used to draw a sample from an arbitrary distribution using a uniformly distributed random number. Some sampling methods are presented in Section 5.2.

Variance, Error and Convergence Rate

The variance of a random variable X, noted V ar(X), is dened as the square value of the standard deviation σ(X) and is computed as follows:

V ar(X) = σ 2 (X) = E[(X -E[X]) 2 ] = E[X 2 ] -E[X] 2 (1.9)
The error of the Monte-Carlo estimator N is dened as:

N = Q N -I (1.10)
It can be demonstrated that:

V ar(Q N ) = 1 N V ar(Q 1 ) ⇔ σ(Q N ) = 1 √ N σ(Q 1 ) (1.11)
which proves that the convergence rate of a Monte-Carlo estimator is in O( √ N ). Thus to divide the variance by two, the number of samples must be multiplied by four. Another way to reduce variance is to reduce the variance of Q N , this is exactly what importance sampling does.

According to [START_REF] Kalos | A Bit of Probability[END_REF], the variance σ 2 N of the estimator Q N can be estimated using the following equation:

σ 2 N = 1 N -1 N i=1 (f (x i ) -Q N ) 2 = N N -1 1 N N i=1 f (x i ) 2 -Q 2 N
(1.12) Thus an estimator of the variance of the estimated mean is given by

V ar(Q N ) ≈ 1 N -1 1 N N i=1 f 2 (x i ) -Q 2 N (1.13)

Example: Estimating the Value of π

We can use Monte-Carlo integration to estimate the value of π. Consider a square which length size is 2 and its inscribed circle C (see Figure 1.7). If we uniformly and randomly create points in the square, the ratio of the number of samples inside C over the total number of samples converges to Area C Areasquare = π 4 . In this case, Monte-Carlo integration is used to estimate the area of the circle C, which is known to be equal to π, with a constant pdf equal to

Importance Sampling

As opposed to uniform sampling, importance sampling tries to maximize the Monte-Carlo estimator by drawing more samples where the value of f is high.

It is a well known variance reduction technique. In fact the perfect sampling would be the one that gives pdf (x) = c * f (x), with c a constant: c = 1 Ω f (x)dx , which leads to:

lim N →∞ Q N = lim N →∞ 1 N N i=1 f (x i ) c * f (x i ) = 1 c = Ω f (x)dx = I, (1.14)
giving the estimator a variance of zero. But, this strategy is only possible when the value of the integrand is known in advance.

Stochastic Ray Tracing

In the previous sections, we have shown that the rendering equation can provide a solution to compute global illumination and we have also explained how Monte-Carlo integration can be used to solve the rendering equation. We now present dierent algorithmic solutions that exist in Computer Graphics to put into practice these mathematical tools, starting with stochastic ray tracing and its derivatives.

To generate the image, the ray tracing algorithm "launches rays", dened as a pair of a 3D starting point and a 3D direction in the 3D scene. For the pinhole camera model, rays start at the camera position and pass through the image plane as shown in Figure 1.8. Then, rays traverse the 3D scene to nd the closest intersection.

Path Tracing

Path tracing (cf. Figure 1.9) is a recursive algorithm based on ray tracing. It launches rays from the camera through the image plane, as described in the previous section, and, each time a ray intersects a surface, launches a new ray from that surface, building a path of light in the 3D scene. A path is thus a sequence of 3D positions that ends when a light is reached or when a stop criterion is attained. We describe these criteria, such as Russian roulette, later on in this document.

Light Tracing

The path tracing algorithm operates in a reverse order compared to what actually happens in the real world. In reality, photons are emitted from light sources and bounce on objects until they get absorbed or reach our eye or a camera sensor. One can render images in the same way by launching rays starting from the light sources and making them bounce on objects in the 3D scene until they reach the camera sensor. This process is call Light Tracing (cf. Figure 1.9). Even though it can render the same images than with recursive ray tracing it is far less ecient. One simple explanation is that the camera sensor represents only a small fraction of the scene, and so, the probability to nd a path that connects a light source to the camera is small. Conversely, light sources represent a larger part of the scene, the path tracing algorithm exploits that property.

16

Arthur Dufay

Explicit Light Source Connection

To ensure that a path is connected to a light source, as shown in Figure 1.8, one can, at each bounce of a ray, make a direct connection. A visibility test needs to be done along this connection to check if light is propagated. This is why such connection is called a shadow ray. This technique is often cited as "next event estimation" or "next event simulation" in the literature, and has been proposed by [START_REF] Kajiya | The rendering equation[END_REF]. By ensuring that the path connects to at least one light source, this technique increases the eciency of the Monte-Carlo estimator.

Several solutions exist to select the light source to connect. The simplest one is to choose deterministically the light source in a round-robin fashion. One can also draw a random number and select uniformly a light source in the set of light sources of the scene. Finally, an ecient solution is to build a cdf over the set of light sources at scene opening, then at each bounce, draw a random number and choose a light to connect according to its potential contribution (i.e., its power) by using the constructed cdf and the associated pdf.

We present a solution, that has been submitted as a patent, to accelerate the computation of shadow rays needed by next event simulation in Chapter 3.

Bidirectionnal Path Tracing

Looking at the two algorithms previously introduced, path tracing and light tracing, we observe that paths starting from the camera sensor (camera paths) may have diculties to reach the light sources, and reversely, paths starting from the light sources (light paths) hardly connect to the camera. Bidirectional 1.4. Stochastic Ray Tracing algorithms, such as Bidirectional Path Tracing (BDPT), try to solve this problem by constructing both camera paths and light paths and connect them at their edges to build more complex paths, as shown in Figure 1.10. Indeed BDPT provides most often a faster convergence rate than path tracing, but it induces a more complex GPU code, and so a slower sample rate. We will see further in detail in Chapter 3, how code divergence and kernel implementation on the GPU is the key to maintain a good eciency. 1. Theoretical Background the specular material BRDF that produces them (cf. Figure 1.11). In fact, as paths are traced from eye to light in path tracing, it is dicult to build a LS + DE path. We present solutions to this problem in Section 1.6, such as the photon mapping algorithm for instance.

As ray tracing requires a lot of computation, another alternative for rendering exists: rasterization that we describe in the next section.

Rasterization

Ray tracing requires a lot of computing power. An alternative algorithm to generate images exists, rasterization and its rendering pipeline (cf. Figure 1.12). Rasterization produces images from 3D models composed of quads or most of the time triangles by projecting them on a 2D plane. It can be executed on a CPU or a GPU, but, thanks to its simplicity and scalability to highly parallel architectures, it has been established as the standard for realtime rendering algorithms on GPU over the years. In fact, rasterization has promoted the use of GPU for real-time graphics.

Figure 1.12: A simplied view of a GPU rasterization pipeline.

Rasterization requires a set of matrices MVP. M represents a model matrix, that transforms a triangle from object coordinate space to world coordinate space This matrix is used to move, scale and rotate objects in the scene as desired. V is a view matrix, that denes the camera position and orientation (i.e., the viewpoint for a particular rendering). Finally, P is a projection 1.5. Rasterization matrix, that denes the view frustum of the camera. Figure 1.13 shows a perspective view frustum.

Figure 1.13: A perspective view frustum.

Fast Removal of Invisible Geometry

To further accelerate rasterization, some techniques are used to discard data that do not contribute to the nal image, we present them here.

Culling and Back-Face Culling

The culling step is done before applying the set of transformation matrices to a triangle. It consists of rejecting all triangles that are totally outside of the viewing frustum. An optional back-face culling operation can also be used. In fact, as 3D scenes are made of 3D objects, which often have a thickness, it is impossible to see a triangle from behind because there is always a front face closer to the camera. Therefore, activating back-face culling can save computation time by discarding all triangles that have a dot product between camera-to-triangle and triangle normal greater than 0 (see Figure 1.14).

Clipping

After transforming a triangle in image space via the transformation matrices, a triangle may fall partially or totally outside of the image. The process of discarding this triangle out of the image is called clipping.

20

Arthur Dufay

Z-Buer

The rasterization process is coupled with the so-called "Z-Buer algorithm". Indeed, when projecting a triangle in the 2D image plane, multiple triangles may fall on the same pixel. To know which triangle has to be stored in that pixel, the Z-Buer algorithm stores the depth of the projected point. Each time a triangle falls into a non empty pixel, its depth is compared to the depth stored. If the new triangle has a smaller depth (i.e., it is closer to the camera), the value in the Z-Buer for this pixel is overwritten. Otherwise, the triangle is discarded.

Rasterization Pipeline and Shading

To this end, we have not discussed what is actually stored in each pixel of the rendered image. We explained how the Z-Buer algorithm works to select the closest visible triangle, but we did not detail the pixel value. To further understand the process of rendering we have rst to introduce the rasterization pipeline.

The basic rendering pipeline shown in Figure 1.12 is composed of two main stages: vertex (resp. fragment) processing handled by the vertex (resp. fragment) shader. These two shaders are two dierent GPU programs, that can be built in the GPU or programmed by the user.

Rasterization

The vertex processing does all the work described in the previous sections to got from a 3D triangle-based model to a set of pixels, also called fragments. It operates on triangle vertices. The fragment shader is invoked for each fragment that actually contains a projected triangle. Its job is to "shade" each pixel of the image. The process of shading is to compute the appearance of a projected triangle in a pixel. To compute shading, a set of variables is passed through the pipeline from the vertex shader to the fragment shader, such as surface normals, 2D or 3D positions, colors, etc.

Shading can be really simple, such as a at shading, or a more advanced one such as the Phong shading (cf. Figure 1.15). The main dierence between these two is that values that passed through the rendering pipeline from the vertex shader to the fragment shader of portions of the triangle between get interpolated in the case of the Phong shading. In at shading only one value for each fragment is given to the fragment shader, typically the value of the rst vertex of the triangle.

Forward vs Deferred Shading

The shading algorithm presented in Section 1.5.2 and shown in Figure 1.12 describes what is called forward shading. In forward shading objects are rasterized and shaded at the same time, one object after another, the Z-Buer taking care of keeping the relevant fragment in the resulting image. Another option is to use a deferred shading pipeline (cf. Figure 1.16). In deferred shading all the objects are rasterized in a rst pass and a set of parameters for each visible triangle is kept in a temporary buer called a G-Buer, the G stands for geometry. Once this is done, a second pass generates an image from the values stored in the G-Buer by shading the fragments. Typical data found in a G-Buer are surface normals, 3D position, material identier.

When the shading is costly to evaluate, the deferred shading algorithm gives better performances than the forward shading one, because it evaluates 22 Arthur Dufay the shading only once per fragment. In comparison, with forward shading, when N triangles project themselves in the same fragment, the shading is evaluated N times. We will see in Section 2.2 how a G-Buer can be used to accelerate path tracing on GPU.

Shading Limitations

The fragment shader can compute simple or more complex appearances but is always restricted to direct illumination or fake indirect illumination, because when it computes the appearance of a pixel it does not have access to the whole scene geometry. We describe in more details direct illumination, indirect illumination and global illumination in Section 1.4. Furthermore, simulating the physics of light and computing a real appearance in a fragment shader would require too much computing resource, and this is not the purpose of a fragment shader.

There are other limitations to rasterization. For instance, it cannot compute real multiple refraction. Some work has been done using textures storing back-faces of objects and nearby geometries (see Wyman [2005a] and Wyman [2005b]). It gives a good approximation and works in real time but it is limited to two interfaces. To compute full refraction, with multiple interfaces, and no approximation on the refracted vectors, the ray tracing algorithm is required. • Finite Element Methods (Radiosity)

• Precomputed Radiance Transfer (PRT)

• Photon Mapping (PM) and its extensions: Progressive Photon Mapping (PPM), Stochastic Progressive Photon Mapping (SPPM)

• Instant Radiosity also known as Many lights methods or Virtual Point Lights (VPLs)

• Monte-Carlo Ray Tracing

• Bidirectional Hybrid Algorithms

In this section we briey describe them to further argument our choice for the path tracing algorithm as a solution for interactive previsualization of VFX. The state-of-the-art report of [START_REF] Ritschel | The state of the art in interactive global illumination[END_REF] is the starting point of our study. They widely cover all the global illumination techniques.

Finite Element Methods

Finite elements methods, also call radiosity, were rst introduced to Computer Graphics by [START_REF] Goral | Modeling the interaction of light between diuse surfaces[END_REF]. To compute global illumination, they rely on a set of geometrical patches that discretize the scene surfaces. Every patch stores a precomputed radiosity value. The main advantage of this method is that it allows fast camera movements. Indeed, once the values of the patches have been precomputed, solving the rendering equation is done easily by fetching the values in the surrounding patches.

More recently, Thiedemann et al. [2011] used voxels to store precomputed irradiance. Even though their method can render global illumination quite fast on the GPU they are limited to two-bounce global illumination.

Despite the quick render time this kind of methods oers, we found that they require too much precomputation and memory. We did not consider them as a good solution in our context.

Precomputed Radiance Transfer (PRT)

Precomputed radiance transfer methods (PRT) were rst introduced by Sloan et al. [2002]. They used spherical harmonics (SHs) to store a transfer func-tion that includes both shading and visibility. As SHs can encode only lowfrequency functions, their method is limited to diuse shading or mid-glossy shading, but it is fast. Indeed, shading computation only requires a dot product.

PRT methods have been extended using other bases to encode them. For instance using wavelets by [START_REF] Ng | All-frequency shadows using non-linear wavelet lighting approximation[END_REF]. A good survey of PRT methods can be found in [START_REF] Ramamoorthi | Precomputation-based rendering[END_REF].

To support highly specular surfaces PRT methods have been extended using spherical gaussians (SGs) by Wang et al. [2009a], and anisotropic spherical gaussians (ASGs) to handle anisotropic BRDFs. Their method provides good result in real time, but are still limited to static scenes.

Finally, Xu et al. [2014] introduced a new method to compute interreections using SGs. Their method works in real time, but is limited to one-bounce interreections, and so, it cannot truly compute global illumination. Due to their limitation on the BRDFs they can handle in some cases, their precomputation step, and their restriction to static scenes, we did not consider these methods as a good solution for our previsualization tool.

Photon Mapping

The Photon Mapping (PM) algorithm was introduced by [START_REF] Jensen | Global illumination using photon maps[END_REF]. As the light tracing algorithm it launches paths from light sources. However, instead of trying to reach the camera sensor, it stores photons at each vertex of the light path in a dedicated kd-tree (the photon map). Then, in a second pass, called gathering pass, the contribution of photons that are in a neighboring area of the shaded surface is accumulated. The main advantage is that this method is more ecient to generate caustic paths LS + DE, as paths start from light sources. The initial algorithm proposed by [START_REF] Jensen | Global illumination using photon maps[END_REF] is biased, but has been improved by [START_REF] Bibliography Hachisuka | Progressive photon mapping[END_REF]. In fact by having a progressive algorithm, that decreases the photon gathering kernel size progressively, [START_REF] Bibliography Hachisuka | Progressive photon mapping[END_REF] ensure that the bias converges to zero, and so make the integrator consistent.

More recently, Stochastic Progressive Photon Mapping (SPPM) introduced by [START_REF] Hachisuka | Stochastic progressive photon mapping[END_REF], also improves PPM by adding the ability of computing radiance over a region instead of only a point. Indeed, the PPM method was restricted to the computation of radiance at a point, rays starting from that point progressively average the contribution of photons. The SPPM algorithm is then capable of computing more complex ray tracing eects such as depth of eld and anti-aliasing.

Despite its robustness, the rst output of the progressive photon mapping (PPM) can have a strange look. In fact they tend to generate some ugly artifacts when the number of gathered photons is too low (see Figure 1.17). Furthermore, the algorithm is not artist friendly in our opinion, parameters Figure 1.17: Artifacts of Photon Mapping when the number of photons is too low. Rendered using Mitsuba, 2500 photons, 119 seconds on an Intel i7-4790K 8 cores CPU. like the gathering radius kernel or the number of photons emitted, are not adjusted easily by a VFX artist. Even though PPM performs really well in LS + DE, it is also known to be slower than path tracing for large outdoor scenes where most of the illumination comes from the sky. For all these reasons we did not choose PM or PPM for our previsualization tool.

Many Lights

The many lights, instant radiosity, or virtual point lights (VPLs) methods were rst introduced by [START_REF] Keller | Instant radiosity[END_REF]. Like the PM algorithm it is a two-pass algorithm that starts by launching rays from light sources. It stores VPLs at each vertex of the light path instead of photons. In a second pass, every VPL is considered as a point light source that emits light uniformly in all directions. It then permits the reutilisation of VPLs for every surface of the scene, decreasing memory consumption compared to the PM solution.

The many lights methods have been improved over the past decades several times. A good overview is given by Dachsbacher et al. [2014a]. It reviews all the techniques that make VPLs scalable. For instance when dealing with a large number of VPLs, one can use Matrix Row Column Sampling introduced by Ha²an et al. [2007]. Their technique uses a matrix where columns store VPLs and rows pixels to shade. By shading only a few surfaces (rows), they can detect the most relevant VPLs (columns). Then, after a clustering step, they can shade all pixels with a reduced number of VPLs.

Other techniques that try to cluster the VPLs exist, such as Lightcuts introduced by Walter et al. [2005] and the more recent Bidirectional Lightcuts by [START_REF] Walter | Bidirectional lightcuts[END_REF]. They both rely on a hierarchical tree structure that organizes VPLs in such a way that only a small number of all VPLs is needed to shade pixels. However, most of these techniques rely on Shadow Maps to compute visibility between surfaces and VPLs, which consumes a lot of memory when dealing with a large number of VPLs. Furthermore, to support specular materials, a lot of VPLs must be shot. In addition to that, the rst iterations of VPL based algorithms tend to generate "splotches" (see Figure 1.18), giving a resulting image that is too dierent from a nal render image in our opinion. This is due Figure 1.18: Artifacts of VPLs often called "splotches" due to a close distance between a VPL and a surface. to a 1 d 2 factor, with d the distance between a VPL and the shaded surface in the computation of the VPL contribution that generates high energy. For these reasons, we did not consider them as a good solution in our context. However, we will see in Chapter 5.5 that VPLs and path tracing can be combined to take benet from both strategies.

Monte-Carlo Ray Tracing

As described in Section 1.4, Monte-Carlo ray tracing methods solve the rendering equation by launching rays in the 3D scene. They can be classied in three types:

• Unidirectional ray tracing: Path Tracing (PT) and Light Tracing (LT)

• Bidirectional path tracing (BDPT)

• Metropolis light transport (MLT)
The unidirectional methods cover path tracing and light tracing. As previously explained, they are easy to implement but do not perform well with complex light paths.

Conclusion

Bidirectional path tracing performs better but is harder to implement efciently on the GPU. In certain cases, for instance outdoor scenes, it is also outperformed by path tracing.

Finally Metropolis light transport (MLT), introduced by Veach and Guibas [1997], is probably the best algorithm to solve complex light paths but in most cases, it has poor performances. It is based on a mutation strategy that, once a light path has been found tries to slightly alter it to nd new light paths that have a high contribution to the image. A good practical introduction to MLT can be found in [START_REF] Cline | A practical introduction to metropolis light transport[END_REF]. Due to its complex heuristic, it is on average outperformed by both PT and BDPT.

They are in our opinion the more versatile methods, they do not require any precomputation, except for the spatial acceleration data structure (cf. Section 4), and are highly parallelizable. The code of a path tracer can be well tailored to t in the GPU (cf. Chapter 3) and is easily adapted to fully use the GPU computation power. For all these reasons we think it is the best algorithm in our context to solve the rendering equation.

Bidirectional Hybrid Algorithms

More recently some hybrid algorithms have been introduced, for instance the Vertex Connection and Merging (VCM) by Georgiev et al. [2012a]. It is a hybrid technique that combines both BDPT and PPM in a nice unbiased algorithm. Although their algorithm can generate complex light paths, it is not straightforward to implement on a GPU. Furthermore, to fully understand their solution, one has to start by implementing a path tracer. In fact, it is a good starting point to learn how to write GPU algorithms properly. For this reason we started with the implementation of a path tracer, with in mind that the next step would be a bidirectional solution like BDPT or a hybrid solution such as VCM.

Conclusion

In this chapter we presented mathematical tools needed to render an image. We also reviewed the dierent Computer Graphics methods to solve the rendering equation and to compute global illumination. Based on what we presented, the path tracing algorithm seems to be a good solution for our VFX previsualization software. It is indeed a versatile, easy to setup and robust algorithm that can fulll our requirements.

The industrial context of this thesis is presented in the next chapter as well as practical Computer Science tools needed to implement our path tracer.

28

Arthur Dufay

Chapter 2

Proposed Path Tracing

Architecture in 3DCast

As previously explained in Section 1.1.1, VFX artists need a software solution to obtain a previsualization of the VFX they are designing. This is especially true at the lighting stage, when light sources, materials, and lighting eects are set up by lighters. This chapter is dedicated to the industrial context of this thesis. We present here a solution that Technicolor provides to artists: the 3DCast platform (cf. Section 2.1). This mixed reality platform is the backbone of all the work that has been done during this thesis. It is able to render complex 3D scenes featuring massive lighting or volumetric rendering.

Unfortunately, the 3DCast platform does not support global illumination, hence the need to extend it with a path tracing solution described in Section 2.2 as well as the dierent features (materials, light sources, ...) that are supported by our rendering engine.

One crucial point for a previsualization software is interactivity. We present in Section 2.2.4, our solution to achieve a better interactivity in our hybrid GPU path tracer. 29 2.1. 3DCast

3DCast

At Technicolor, the VFX Interactive Synthesis Team, develops and exploits 3DCast, a mixed reality framework that enables real-time visualization of 3D virtual worlds on networked devices such as mobile phones, tablets, desktop or laptop PCs. It allows real-time animation and rendering of complex 3D virtual worlds (see Figure 2.1). The topics covered by the research and engineering teams using the 3DCast platform include: progressive meshes, procedural models, human animation, 3D interfaces, facial expressions, volumetric rendering (Figure 2.2) and massive dynamic lighting (Figure 2.3), among others. Within 3DCast, worlds can be fetched from local drives or streamed through heterogeneous networks, including LAN, ADSL, Wi. Furthermore, the animation, interaction and rendering can be performed on PC clusters, workstations, or smart-phones (i.e., heterogeneous clients). For performance reasons, only a subset of the virtual world visible from the current viewpoint may be transmitted and visualized. In this case, when navigating through the virtual world, streaming algorithms may anticipate which elements of the environment need 30 Arthur Dufay to be streamed to the client. Hence these algorithms combine low bandwidth consumption with high reactivity. 3DCast is a modular high-performance virtual reality platform based on an extended X3D scene graph as well as a set of plug-ins: OpenGL scene graph implementation, streaming of large terrains, video eects, etc.

Based on this description, the platform is made of three abstraction layers:

• The system layer provides generic encapsulation of system dependent calls. This ensures the multi-OS interoperability of the platform.

• The application layer provides generic containers and tools for application creation. It also provides some generic built-in components:

The distributed meta scene graph and Internet protocol allow for 3D world distribution and synchronization.

The generic scene graph renderer allows for heterogeneous 3D world adaptive rendering (see Figure 2.4).

• The components layer conceptually regroups all the components (builtin and "external" ones). Components can be added at compilation time or at run-time through the use of a plug-in system provided by the application layer. Thus, developers can quickly design their own applications for the platform. For instance, the current component set includes: VRML/X3D scene graph, a plug-in for large landscapes, several video eects, high denition textures management, etc. These built-in components provide generic data structures and algorithms for the generated 3DCast application. They also ensure interoperability of the "external" components.

The 3DCast platform is delivered to Technicolor clients and associated companies, like The Moving Picture Company (MPC), as both a standalone software and an Autodesk [2017] Maya plugin. The 3DCast renderer plugin for Maya really helps lighters by giving them a better render than the standard Maya Viewport. It provides them a better feedback, for instance using our contact visualization system, see Figure 2.5, that helps them to place objects in 3DScene. This tool was published by [START_REF] Marvie | Contact Visualization[END_REF]. Even though 3DCast can render complex lighting eects, it cannot compute global illumination. This is a well-demanded feature by lighters at the previsualization stage. In fact, giving them the ability to previsualize what the nal render can look like in a few seconds or minutes instead of a few hours using a nal render algorithm is a tremendous asset. The main challenge of my PhD was to enhance 3DCast with a global illumination solution. We made the choice of focusing on a path tracing solution because it is, in our opinion, the most versatile global illumination algorithm that simulates light transport without the bias that can be introduced by rasterization techniques. It is also the solution that would give the closest result to the nal render image, because most of the nal render production renderer uses path tracing. Table 2.1 shows the dierent algorithms used by some well-known production renderers. Path tracing is still today the "gold" standard algorithm used in production even though it does not perform well when rendering complex light paths such as caustics as explained in Section 1. in Section 1.6.6, start to appear in production, for instance at Pixar [2017] in their Renderman renderer. Since we chose to focus our work on path tracing, our main goal is to make it faster using the GPU, and doing so applicable on a lighter personal computer for previsualization of VFX.

Path Tracing in 3DCast

We describe in this section all the path tracing features that are, in our opinion, mandatory for a previsualization path tracing engine and thus, were implemented in 3DCast during this PhD. Obviously, this is still a work in progress and some features are lacking, they are documented as future work in Section 2.3.

GPGPU

The rst task to develop a GPU path tracer was to support General-purpose processing on graphics processing units (GPGPU) in 3DCast. It permits to launch any kind of computation on the GPU using computation kernels. We will not enter in detail of the GPU compute capabilities in this section, since Chapter 3 is dedicated to that. Several APIs for GPGPU are available:

• Direct3D Compute Shaders • NVIDIA CUDA • OpenGL Compute Shaders • OpenCL
As we did not want to be restricted to Windows users, so we rejected the Direct3D option. At the the time we started working, OpenCL was slower than OpenGL Compute Shaders or NVIDIA CUDA. It was mostly due to memory transfer between OpenGL context and OpenCL context, so we rejected it also. That left us with NVIDIA CUDA and OpenGL Compute Shaders. As 3DCast was already using OpenGL for rendering, and NVIDIA CUDA being limited to NVIDIA graphic cards we chose to implement our prototypes using OpenGL Compute Shaders. With hindsight, NVIDIA CUDA would have been a better choice. In fact, being limited to NVIDIA graphic cards is not a big issue. Furthermore, NVIDIA provides some very powerful performance analysis tools with the CUDA framework. To this end, the proling of OpenGL Compute Shaders is still at its early stages.

Materials

Four BRDFs are implemented in the 3DCast Path Tracer. We think that these four models are sucient to represent a wide variety of materials encountered in 3D scenes.

Normalized Phong BRDF The normalized Phong BRDF, introduced by Lafortune and Willems [1994], whose formula is:

brdf P hong = k d π × C dif f + k s (e + 1) 2π × C spec × (r • ω o ) e ω i • n (2.1)
with n the normal at the surface, ω o and ω i respectively the outgoing and incoming vector and r the reection of the vector ω i over n. e controls the shininess of the material, a small value (under 10) gives a rough appearance to the material, whereas a large value (over 1000) gives a highly specular appearance. k d and k s are oating point values in [0, 1], with these Perfect Mirror A perfect mirror BRDF is also implemented, it can be represented as a Dirac, that reects light in a unique direction, the reected direction of light.

brdf mirror = k s × C spec × δ(r -ω o ) = k s × C spec if r = ω o 0 otherwise.
Refractive materials Refractive materials are also handled, using an approximation introduced by [START_REF] Schlick | An inexpensive brdf model for physicallybased rendering[END_REF]:

R schlick = R 0 + (1 -R 0 )(1 -ω i • n) 5 R 0 = n 1 -n 2 n 1 + n 2 2 (2.2)
n 1 and n 2 are the indices of refraction of the two media at the interface. R schlick approximates the Fresnel term: it gives the ratio of reected/refracted light.

Path Tracing in 3DCast

Cook Torrance A more advanced BRDF is also implemented, the Cook Torrance BRDF introduced by Cook and Torrance [1982]. We compute it using the following formula:

brdf CookT orrance = k s × C spec × F × G × D 4(n • ω o )(n • ω i ) (2.3)
We use it with a Beckmann distribution, where the D term corresponds to the microfacet distribution, F is the Fresnel term and G is the geometrical attenuation. We use a [START_REF] Schlick | An inexpensive brdf model for physicallybased rendering[END_REF] approximation for the F term. D, F and G are computed as follows:

D = e (n•h) 2 -1 m 2 (n•h) 2 πm 2 (n • h) 4 G = min 1, 2(n • h)(n • ω o ) (ω o • h) , 2(n • h)(n • ω i ) (ω o • h) F = F 0 + (1 -F 0 )(1 -(ω o • h)) 5
(2.4) h is the halfway vector and computed as:

h = ωo+ω i |ωo+ω i | .
The m parameter ranges in [0, 1] and controls the roughness of the material. F 0 is the material response at normal incidence, it can be computed from the refractive index of the material µ using the following formula:

F 0 = 1-µ 1+µ 2 .
Importance sampling Following the equations given in Dutré et al. [2001], importance sampling is implemented for the Phong BRDF. Samples for a diffuse lobe are generated using the following cosine weighted distribution:

x = cos(2πr 1 ) 1 -r 2 2 y = sin(2πr 1 ) 1 -r 2 2 z = r 2 (2.5)
where r 1 and r 2 are random numbers in the range [0, 1]. The pdf of such random direction is pdf (θ) = 1 2π . For a glossy lobe, we use a sampling method proportional to the power exponent of the BRDF:

x = cos(2πr 1 ) 1 -r 2 e+1 2 y = sin(2πr 1 ) 1 -r 2 e+1 2 z = r 1 e+1 2 (2.6)
where e is the Phong exponent, the pdf of such random direction is pdf (θ) = They allow the setup of dierent lighting eects. We present them here. As for material colors (C dif f and C spec ) previously dened, light source colors are also dened as RGB color encoded with 8 bits per color channel, and so, do not feature spectral rendering.

• Point Light: it is dened by its position, its color (in RGB space) and its intensity. It represents a single point in 3D space and thus has no equivalent in the real world. It emits the same amount of light in all directions. Its contribution is evaluated as follows:

L o (x, ω o ) = LightColor * LightP ower[W ] * f r (ω i , ω o )[sr -1 ] ω i 2 [m 2 ] (2.7)
where ω i 2 is the distance between the point light and the shaded point x.

• Spot Light: in our implementation, a spot light is a point light with one restriction. It does not emit light in all directions but only in a cone (cf. Figure 2.8). Spot lights can also be implemented with a decay factor over the cone, having rays farther from the main direction of the spot light emitting less light. We did not take into account this decay factor in our implementation.

Figure 2.8: The Mitsuba cap model from [START_REF] Jakob | Mitsuba renderer[END_REF]. Rendered with a copper material, standing on a plane and lit by a spot light with no decay factor and a uniform blueish environment light.

• Environment Light: An environment light has no position in the 3D scene: its purpose is to illuminate the entire scene. It can be represented as a sphere englobing the whole scene and positionned at innity (cf. Figure 2.9).

Figure 2.9: The mitsuba cap model, from [START_REF] Jakob | Mitsuba renderer[END_REF], with a glossy material, standing on a plane and lit by an environment light.

• Area Light: An area light (cf. Figure 2.10) provides a closer representation of the real world, having a surface and providing smoother and more realistic illumination. In our implementation, we consider area lights as polygonal diuse emitters, limited to rectangular shapes. It emits the same amount of light in all directions that have a positive scalar product with the normal of the area light and nothing in the others (i.e., it emits light in front of its geometry and nothing behind it). Compared to point light or spot light sources they have the advantage of casting nicer soft shadows and give a more realistic look. The drawback is that they need to be sampled: multiple rays have to be sent to compute their contribution. To solve that problem analytic solutions exist, see [START_REF] Lecocq | Accurate analytic approximations for real-time specular area lighting[END_REF] or [START_REF] Heitz | Real-time polygonal-light shading with linearly transformed cosines[END_REF], but none of them take into account visibility and occlusions.

Russian Roulette

When using Monte-Carlo path tracing, a lot of computation time can be spent on launching rays that contribute faintly to the image. In fact, looking at the innite sum of the rendering equation (cf. Equation 1.3), we see that as the path length grows, we add samples with less and less energy because energy gets absorbed along the path by BRDFs.

To solve that problem three potential solutions exist. The rst one is to limit the path length by a xed threshold, this ts well to the GPU but some bias is introduced, due to energy lost by paths that would have been longer if they were not stopped by the threshold.

Another solution is to use Russian roulette, a technique that was introduced 38 Arthur Dufay by [START_REF] Veach | Robust Monte-Carlo Methods for Light Transport Simulation[END_REF]. Even though it may increase variance in some cases, it has the benet of increasing the eciency of the Monte-Carlo integrator by keeping a constant contribution of samples along a path. This is done by adding a probability p for each path to be stopped at each bounce. We x p to the inverse of the absorption factor of the BRDF. This is a common choice that permits to simplify the computation of the Monte-Carlo integrator. For instance, in the case of a diuse BRDF with an absorption value of 1 -k d , we set p = 1 k d . In the Monte-Carlo integrator, p and k d of the BRDF get canceled out, leaving us with a constant energy along the path.

Still, using Russian roulette, paths could be very long. This is why we opted for another solution. We combined Russian roulette with a maximum xed length of path. By setting this maximum value very high we minimize the bias and still save computation time.

Tone Mapping

All our rendered images are stored using 32 bits per channel corresponding to the OpenGL format RGBA32F. This is mandatory since our rendering algorithm computes images that have unbounded oating values. For instance, a point light source with a power of 1000 watts that illuminates a surface with a diuse BRDF of k d = 0.5 can lead to a pixel value of 500. Furthermore, we compute our images iteratively. Indeed, the results of the Monte-Carlo integrator are added at each frame, as described in Section 2.11. Hence the need to store the results in a high precision buer as we do not want to loose information. However, a typical computer screen has a limited luminance range. It displays images with only 8 bits per channel. To display our images, a tone mapping operator converts our high dynamic range images (HDR) to low dynamic range images (LDR). We used the following linear tone mapping operator, that applies a gamma correction and rescales luminance: When dealing with path tracing on GPU, each time a ray is launched to nd an intersection, the entire 3D scene needs to be accessible to avoid starvation of threads waiting for geometry to test. It is worsened by slow data transfer between GPU memory and CPU memory. Fortunately, NVIDIA provides an extension to address the whole GPU memory with data pointers. The GL_NV_shader_buer_load, see [Brown et al., 2010], provides a mechanism to organize pointers and fetch any data on the GPU from a shader. Using this, we can access all the geometry of the scene as well as its corresponding spatial acceleration data structures, described in Chapter 4, each time a ray is launched by a compute shader. Obviously, this works as long as the 3D scene ts in the GPU memory. We did not develop any out of core solution during this thesis.

Rasterization as Primary Ray

To further accelerate computations in our path tracer we developed a hybrid path tracing pipeline, using a G-Buer from a rasterization pass as the primary rays (i.e., rays starting from the camera). This pipeline is described in Figure 2.11. Using the G-Buer we can start the path of our path tracer at the surface stored in the G-Buer. To keep an interactive framerate the path tracer computes paths in an iterative manner, launching only one path per pixel per frame. At each frame, the path tracer launches rays until the path is terminated, then its contribution is added to the image as shown in Figure 2.11.

Figure 2.11: Our hybrid path tracing pipeline, using a G-Buer from a rasterization pass to replace primary rays. Top red square shows the rst pass that computes the G-Buer using a rasterizer. Middle red square shows a pass of path tracing. Bottom blue square shows the full algorithm composed of one pass of A and several passes of B.

However, this method has some drawbacks, for instance, we can not simulate ray traced depth of eld. To do that we would have to replace the G-Buer by a full path tracer, launching rays starting from the camera, and therefore simulating the camera lens behavior. Using this G-Buer technique also introduces some spatial bias.

Quad-Tree Pixel Sampling

When trying to compute global illumination in a complex scene, even launching only one path or even just one ray per pixel can take several seconds on the GPU. Thus, to reach an interactive framerate, we can compute only a subset of the image at each frame. In fact as the launching of computation kernels on the GPU blocks any interaction that the user might have with the computer, the computing time must be as short as possible. This is a really important point in the context of a previsualization tool. It must give feedback to the user as fast as possible without interrupting its work.

To reach this interactive framerate we chose to divide our image into square tiles of 256 pixels (i.e., tiles of 16 pixels per side). At each of these 256 subframes, we compute only one pixel in each tile. Note that this could have been extended to tiles of adjustable size.

Figure 2.12: Our 256 iterations pattern for lling a path tracing tile of 16×16 pixels.

The order in which we choose the pixel to compute at each sub-frame follows a hierarchical structure similar to a quad-tree. For the four rst subframes, we select the four left-top most pixels in each sub-tile of 4×4. The 12 next frames: frames 5-16, compute the 12 pixels in each sub-tile of size 2×2, and so on. This order is presented in Figure 2.12.

By doing that, we can splat the value of already computed pixels in the tile to non-computed pixels really fast We use a simple loop that starts at the higher resolution of the quad-tree and goes to a lower resolution until it nds a computed pixel, see Algorithm 2.1. Note that this algorithm splats only values for the rst frame, from frame 2 all the pixels in the image will be computed at least once, so splatting will not be needed anymore.

Algorithm 2.1: Splatting algorithm for the rst frame of our path tracer with a tile of 256 pixels. // otherwise go to next resolution factorSize = factorSize * 2; texelPos2 = ivec2 ( int ( texelPos . x / factorSize ) , int ( texelPos . y / factorSize ) ) ; texelPos2 = ivec2 ( int ( texelPos2 . x * factorSize ) , int ( texelPos2 . y * factorSize ) ) ; } This splatting technique also helps us to have a quicker feedback on indirect lighting when moving the camera or editing the scene, see Figure 2.13.

Even though this multiple frame computation improves interactivity it has a major drawback. The computation time required to compute a full image at 1 sample per pixel (spp) takes more time than a full image computation at 1 spp in one frame. We think that this is mostly due to OpenGL driver overhead and GPU cache misses. Indeed, when we split the computation of an image into multiple sub-frames, pixels that are close to each other in image space will not be computed at the same time with our quad-tree repartition. These pixels tend to generate similar rays in 3D space and so have some GPU pixels on an NVIDIA GTX970, using three dierent methods. The one-pass method computes the image in one GPU kernel launch. The 16-pass blocks methods launches 16 passes of GPU kernels on 1 16 of the entire image. And nally our 256-pass quad-tree described in Section 2.2.4 and Figure 2.12. The "Time Ratio" column shows the performance factor over the fastest method (the one-pass method). Note how the 16-pass method as well as the 256-pass quad-tree method decrease performance. The six test scenes are presented in Figure 2.14. cache coherence. The full image computation method is then more performant than our multiple frame computation, it computes neighboring pixels in the same GPU warp, at the same time. We conducted a benchmark to measure the overhead of our method, results are shown in Table 2.2. We added in this benchmark another method of rst frame computation: the 16-Pass blocks method, which computes the image in 16 passes on the GPU by operating on each pass on a sub-frame of 1 16 of the entire image. This benchmark also shows the increase interactivity we obtain using our quad-tree method (see the FPS column).

Conclusion

We have presented in this chapter the applicative context of this thesis as well as the features we have implemented in our path tracer. This basis helped us to better understand the challenge of implementing a production renderer.

Obviously, there are still some features lacking, as future work for the materials. We would like to implement any material that can be used to simulate skin behavior, for instance using a bidirectional scattering-surface As explained in Section 2.2.4 we can access the whole GPU memory but we do not provide any mechanism to handle scenes that do not t on GPU memory. Thus some future work would be to address out of core techniques.

Finally, in Section 2.2.4, we provide a solution to obtain a better interactivity in our hybrid rendering pipeline.

In the next chapters, we will further explain where the challenges were to implement this path tracer. For instance by analyzing how to properly write GPU kernels, and take care of potential cache misses induced by kernel code (cf. Chapter 3). We will also see how to accelerate algorithmically ray queries using spatial acceleration data structures in Chapter 4. Finally, we will further increase performance of our path tracer by leveraging GPU computing power (cf. Chapter 5) with our new decorrelation technique.

For a long time, CPU has been the main central computing unit on a computer. However, since the last decade, GPUs and their huge computing power (cf. Figure 3.1) have become more and more used as computing units. For instance the latest NVIDIA GPU for the consumer market, the GTX 1080, has a processing power of 7967 GFLOPS (Giga Floating point operation per second) at single precision and 2560 cores, whereas the latest Intel CPU for mass market, the i7-6900K has 8 cores (16 threads supported with hyper threading) and a computing power of 819 GFLOPS at single precision. The GFLOPS for a CPU can be calculated using the following formula:

GF LOP S CP U = # cores × core f requency (GHz) × OPC (3.1)
where OPC is the number of operations per clock cycle. For an Intel i7-6900K the OPC rises to 16 at double precision and 32 at single precision. For a GPU, the formula is a little bit dierent. GPUs have the ability to compute a mul-add, the combination of a multiplication and an addition, for instance, mul hardware. The GFLOPS on a GPU can be calculated using the following formulas:

a = b × c + d,
OP C GP U = (# mul-add units × 2 + # mul units ) GF LOP S GP U = # cores × # SIM D units × OP C GP U × core f requency (GHz) (3.2)
It is not surprising that highly parallelizable algorithms are now set up on GPUs, but, to prot fully of the GPU computing power, one needs to take great care when writing computing kernels. GPUs have a very peculiar architecture that must be taken into consideration. A good introduction to GPGPU and GPU architecture can be found in [START_REF] Kirk | Programming Massively Parallel Processors: A Hands-on Approach[END_REF].

GPU Architecture

We now introduce in detail the GPU architecture and its components. We base our presentation on the NVIDIA Kepler architecture described in [NVIDIA, 2012a] and [NVIDIA, 2012b]. Note that more recent architectures have been introduced since: the Pascal and Maxwell architectures. However the Kepler architecture remains sucient to understand the main principles of the GPU.

GPU Cores Hierarchical Structure

Hardware Structure

An NVIDIA GPU has a hierarchical structure (cf. Figures 3.2 and 3.3). CUDA cores are grouped into warps. Up to this day on any NVIDIA card, a warp is always made of 32 threads. Warps are then grouped together to ll a Streaming Multiprocessor (SM or SMX, X stands for next generation, as opposed to the older Fermi's architecture SM). Each Kepler SMX contains exactly 192 cores (i.e., 6 warps). This is true for every NVIDIA Kepler GPU. To vary the number of cores available on a card, NVIDIA varies the number of SMs. SMs are then grouped into Graphics Processing Clusters (GPCs). The highest level of the hierarchy, that grouped all the GPCs is the grid. Each level of the hierarchy has its own dedicated memory or execution unit. For instance, GPCs have a dedicated raster engine to execute core graphics functions. .

Logical Structure

In addition to this hardware structure, there is a logical structure. When one wants to launch a kernel computation on the GPU, one has to specify both the number of threads and their repartition on the GPU. This is done by specifying the number of blocks that group threads together. Thread blocks are also called workgroups. A memo of the terminology used by CUDA is given in Table 3.1. GPGPU ensures that threads inside a block will be executed together. Thanks to that, some synchronization barriers can be used inside a group. On the contrary, there is no synchronization scheme available for the whole grid, and no assumption can be made on the order in which groups will be executed.

Blocks also allow the use of shared memory, for fast data access and reutilization of data across several threads. After the number of blocks and their size have been specied, the GPU will distribute the blocks on the SMs. Note that all the blocks have the same size. Several blocks can be assigned to one SM as long as the resources available on that SM satisfy the resources needed by the block.

Thread

Lightweight process, to be executed in parallel

GPU Memory Layout

As shown in Figure 3.4, the memory layout inside a GPU is also hierarchical. It goes from fastest to slowest in this order (numbers are given for the Kepler GK104 and GK110): • L2 cache 1536KB on Kepler GK110, 512 KB on Kepler GK104, medium latency (100-300 cycles), high bandwidth ≈ 750 GB/s

• GRAM (also called DRAM) 3GB-6GB on Kepler GK110 high latency (400-800 cycles), low bandwidth ≈ 250GB/s

It starts with the lowest layer the registers. Registers are used to store automatic scalar (non-array) variables. Each thread has a private copy of these scalar variables. Array variables are stored in global memory (GRAM or DRAM), which make them slower to access. In some cases, when using xed size array, the compiler may decide to store an array into registers. The L1 cache also stores the shared memory for the workgroup. On top of that, an L2 cache stores temporarily accessed variables for the grid. Finally, GRAM stores all the data that may be accessed in a kernel code. Even though GPUs have a tremendous computing power, they have several limitations (see below) that must be taken into account.

Memory Access Bottleneck

Memory access is a major bottleneck when using GPGPU. The performance of a kernel execution directly depends on the presence or not of data in the cache. In fact, registers and L1 cache are really fast compared to L2 cache or the even slower GRAM. When a kernel tries to access some data that is not present in the cache, a long fetching operation moves this data from L2 or GRAM to L1 cache. As threads work in parallel inside a warp, increasing the number of memory fetches directly impacts performance. An algorithm 3.3. GPU Limitations that makes the threads work on the same data performs better than one that accesses totally random data.

Register Size Limitation

Another limitation that appears in GPGPU is the number of registers. When there are not enough registers available to fulll kernel requirements, the SMX reduces the number of blocks running in parallel.

Take for example the Kepler GK104. It can have at maximum 2048 threads per SMX, and 1024 threads per block. The register size for its SMX is 65536 (65K) of 32-bit registers, and the number of registers available per thread (if the SMX is fully utilized) is 65536 2048 = 32. If one uses computing blocks of 1024 threads, with 32 registers per thread, the number of threads that can be executed concurrently will be 2048, (i.e., 2 blocks of 1024 threads). Now, with the same thread repartition, if the registers used per thread is increased to 33, there would not be enough registers for 2048 threads to work in parallel on an SMX. The number of concurrent threads will be reduced. Since the reduction in the number of concurrent threads running on a SMX is done at block granularity, in this example, the number of concurrent threads will be reduced from two blocks per SMX to one block per SMX. Only 1024 threads will run in parallel, leading to a 50% of SMX utilization. In this example to increase SMX utilization, the block size must be reduced.

Great care must be taken on the register utilization and the repartition of threads on the grid. This high variation on SMX occupancy shown in the last example happens at some threshold on the number of registers. A good tool to compute the best thread repartition is provided by NVIDIA, the CUDA Occupancy Calculator (cf. [NVIDIA, 2017a]).

Kernel Branching

The Single Instruction Multiple Threads (SIMT) nature of the GPU also imposes some constraints. The GPU executes the same instruction for all the threads. It means that, when using an if-else statement, if at least one thread needs to execute the if-then part, all the threads will execute it with it. Obviously, the same rule applies for the else part. In the case of a for-loop, it is even worse, since all the threads run the loop as long as one thread still needs to run it. This particularity can lead to poor performance. Kernel code with a lot of branching can lead to long execution times. To counter this eect, several solutions exist. One can, for instance, reduce branching by rewriting the code. Another alternative is to split the kernel code (at the branching point) in several sub-kernels. Then the thread pool is split on the CPU to execute the corresponding sub-kernels code on the GPU. The path tracing algorithm may look like a simple algorithm as presented in Algorithm 1, but its ecient implementation on the GPU can be tedious. Several variants of path tracing are presented and well studied in [START_REF] Davidovi£ | Progressive light transport simulation on the gpu: Survey and improvements[END_REF]. In this section, we present some of them and explain from where performance dierences between them might come from. We start by explaining path regeneration principle. Then, we present the two implementations of our path tracer and give a benchmark to compare them.

Algorithm 1 Pseudo code of Path Tracing using a GBuer for primary rays.

1: for all pixels in Image do 2:

for i=0 to nbSamples do

Path Regeneration

When dealing with path tracing on the GPU, one has to face the sparse warp problem: not all threads in a warp are active. This problem comes from both the Russian roulette and the possibility that a path may terminate when it leaves the scene. In both cases, there is a divergence in the length of the paths in the warp. With each thread having its own path to process, the warp occupancy drops as the path length grows. [START_REF] Novák | Path regeneration for interactive path tracing[END_REF] introduce the path regeneration technique, which consists of reassigning a new path to the terminated threads. This strategy maintains a 100% warp occupancy. However, as stated by [START_REF] Van Antwerpen | Improving simd eciency for parallel monte carlo light transport on the gpu[END_REF], this strategy has two main drawbacks. First, during the regeneration phase, there is a code divergence: only threads that need a new path need the regeneration phase. Secondly, it leads to a divergence in the bounce state of the threads inside a warp, i.e., inside a warp, a thread might be computing the third bounce of path while another might be on the rst bounce. This bounce state divergence breaks the primary ray coherence and results in a performance drop. This is even true for secondary rays (paths at their rst bounce) that still have some potential coherence when leaving a at surface, and especially using our decorrelation technique presented in Section 5.4.

To counter this eect, van Antwerpen [2011] proposes a stream compaction in addition to the regeneration. It consists of removing all terminated paths from the stream of threads. The stream of threads is then compacted and the removed samples use the regeneration and are placed at the end of the stream. By doing so, the regenerated samples are executed in the same warp, and SIMD eciency remains high.

First Implementation -Single Kernel Path Tracing

We now present here our rst implementation of our GPU path tracer. All the code is written in one GPU kernel, presented in Algorithm 2. In this version we did not use regeneration, nor stream compaction. This implementation was a good starting point for our GPU path tracer. It is quite similar to the naivePTsk by [START_REF] Davidovi£ | Progressive light transport simulation on the gpu: Survey and improvements[END_REF].

Multiple Kernels

Even though the previous algorithm (cf. Algorithm 2) using a single kernel was already giving good performance, we wanted a more exible pipeline for our GPU path tracer. Indeed, according to Laine et al. [2013], a path tracer implemented on the GPU as one megakernel, is harmful to performance. To solve that problem, they propose a wavefront approach. Using their solution the GPU path tracer is then cut in multiple kernels. They also sort paths on the dierent materials they encounter, and execute a dedicated GPU kernel for each material. This method improves performance when dealing with complex materials, such as car paint which is composed of several layers, and makes it hard to evaluate.

56

Arthur Dufay

Kernel Implementation of Path Tracing on GPU

Algorithm 2 Path Tracing using a single kernel on the GPU. The grid dimensions match the image dimensions (i.e., there is a 1-to-1 mapping between pixels and threads). The kernel is invoked N times to compute N paths. 1: result = vec3(0,0,0) 2: energy = vec3(1,1,1) 3: //setup secondary ray from GBuer 4: //a surface contains the starting point, normal and material 5: fetchGBuer(pathNum, threadID, surface) 6: for bounce = 0 to maxNumberOfBounces do Following that philosophy, we decided to implement a multiple kernel version (see Algorithm 3). It is composed of three GPU kernels.

• The sampling kernel (A)

• The ray tracing kernel (B)

• The accumulation kernel (C) Kernel A, depending on the BRDF of the current ray starting surface, uses importance sampling to generate a new ray direction. It also features the path regeneration, in case the path was terminated by a previous call of Kernel C. Kernel B, launches a ray to nd the next closest intersection. If needed, a shadow ray is shot to account for next event estimation (see Section 1.4.3).

Kernel C is responsible for the Russian roulette termination criterion, as well as storing the path in the rendered image. This implementation scheme is quite similar to the RegenerationPTmk described in [START_REF] Davidovi£ | Progressive light transport simulation on the gpu: Survey and improvements[END_REF]. It has several advantages for our previsualization tool. First, if the ray tracing step in Kernel B takes too much time to compute, due to too complex geometry, we can easily do it in several passes to prevent the GPU from blocking user interaction.

We can also follow the Laine et al. [2013]'s idea, and replace Kernel A by several dedicated kernels based on the material to compute more complex BRDFs and importance sampling techniques.

In comparison to the single kernel implementation described in Section 3.4.2, this version gives us a faster framerate. Indeed, the pipeline computes only a single bounce of the path at each frame, giving a faster feedback.

We also avoid several pitfalls of the GPU using this implementation scheme. First, the number of registers used is smaller (cf. Section 3.3.2), and the number of instructions per kernel is also reduced. Secondly, by using dierent Kernel A based on the material type we can reduce the kernel branching (cf. Section 3.3.3).

58

Arthur Dufay

Kernel Implementation of Path Tracing on GPU

Algorithm 3 Path Tracing using three kernels on the GPU. The grid dimensions match the image dimensions (i.e., there is a 1-to-1 mapping between pixels and threads). The kernel is invoked N times to compute N paths.

Benchmark

We present here (cf. Tables 3.2, 3.3, 3.4 and3.5) a benchmark on the two implementations of the GPU path tracer we presented in the last sections. All these benchmarks were done with a 1280 × 720 pixels viewport. A maximum of 7 bounces was used, except for the San Miguel scene in Tables 3.4 and3 As expected the multiple kernel option is protable. However, its interest is limited by several factors. First, on more recent GPU architectures (Maxwell) the performance gain is smaller than on older ones (Kepler and Fermi). This is due to an augmentation of cache performance, both in size and speed, on recent architectures by NVIDIA. Secondly, on bigger GPUs (more computing kernels) like the GTX TitanX, the performance gain is noticeable only for the largest scene (San Miguel). However, in these test scenes, only a basic Phong shading was tested. We truly believe that we would have a much more protable gain using the multiple kernel option with more complicated materials.

We also think that Kernel B can be optimized, by ray sorting or another option. In that sense, in the next section, we present a patent we submitted to further increase shadow ray queries.

Reverse Shadow Ray

3.6

Reverse Shadow Ray

Technical Problem Solved by the Invention

As presented in Section 1.4.3, to increase path tracing eciency a direct connection to a light source is done at each bounce of the path. However, as the light source is randomly chosen, this leads to a potential shadow-ray divergence. Indeed, GPU threads working in the same warps are gonna try to connect their respective path to dierent light sources.

Our solution helps to reduce computing time in the light sampling stage in path tracing engine. It is based on a reversed shadow ray technique, and helps to get a better coherency between shadow rays computed in the same GPU warp.

Proposed Solution

Instead of launching light visibility rays from the surface to the light, we launch it from the light to the surface. Since those rays are just visibility rays, they only need to return if there is an intersection between two points in 3D space, the light source and the surface being lit, so the orientation they took does not matter.

Reverse Shadow Rays

By doing intersection query in reverse order we improve ray coherency because rays computed in parallel on the GPU will start from roughly the same point in 3D space. In other words the possible ray start position from a light source is greatly reduced compared to a classic light ray query when rays can start from any surface in the 3D scene.

Clustered Shadow Rays

To improve performance we cluster shadow rays by light source origin. By doing that, all the threads in a GPU warp will compute shadow rays starting from the same light source. This clustered solution can be done using several methods such as a fast GPU reduction over a buer containing ray queries. An example is given by van Antwerpen [2011] (referred to stream compaction).

Clustering Algorithm

We consider a 2D grid of threads mapped to 2D grid pixels of the image. We consider a 2D buer as the light sampling buer, also mapped on the computing grid and the image grid. Each thread in the computing grid will be responsible of computing a path for a pixel. At each bounce along the path, the 62 Arthur Dufay

Kernel Implementation of Path Tracing on GPU

next event estimation will randomly choose a light and store its index in the light sampling buer. One can reduce the light sampling buer into N small buers, N being the number of light sources potentially chosen by the path tracer, and then compute reverse ray queries on those N buers to solve the visibility between sampled light sources and surfaces. This clustered method is described in Figures 3.5 and 3.6. Then, it is split in smaller buers to regroup shadow rays per light source.

Advantages of the Method

This method is compatible with any ray tracing system or deferred rendering engine. It is easy to implement and gets faster path tracing without any change to the sampling strategy.

Conclusion

Conclusion

We have presented in this chapter two dierent implementations of a GPU path tracer and review their respective performance. Even though the constant increase in both speed and size of the GPU memory cache limits the potential gain in performance of the multiple kernel implementation, we believe that it is still a good option since the increasing complexity of 3D scenes will counterbalance it. Once again, a more complex surface material than the one we used in our test scenes would have exposed a more protable performance gain. As future work, we would like to implement the stream compaction technique of van Antwerpen [2011]. On a wider eld of research, we would also like to investigate on automatic procedures to implement kernels according to the GPU used.

We also presented our patent Fast visibility test using reversed shadow rays, designed to increase the eciency of shadow ray queries in ray tracing applications. Unfortunately, in our industrial context, it did not result in any ray tracing performance improvement. This is due to a too important platform overhead and an impossibility for us to measure precisely where the caveats of our implementation were. We still believe that it can be a protable option for simpler platforms, and even more using CUDA instead of OpenGL Compute Shaders to nely adjust kernel implementation.

We have demonstrated in this chapter that the consideration of thread repartition, as well as their cache fetching operations, is directly aecting ray tracing performance. In the next chapters we will see how to further increase performance, algorithmically by using dedicated data structures (cf. Chapter 4), and on focusing on sampling strategies (cf. Chapter 5).

Arthur Dufay

Chapter 4

Analysis for an Adequate Spatial Acceleration Data Structure

Since path tracing requires a huge amount of rays to be launched in the 3D scene and thus many ray-mesh intersection tests, it is "mandatory" to organize the scene 3D geometric data in an adequate data structure prior to the rendering stage. In Computer Graphics, many structures known as spatial acceleration data structures (SADS) have been introduced over the past decades. These data structures help to reduce the number of ray-triangle intersection tests by partitioning triangles into subsets of the 3D scene. Theoretically, they can reduce the complexity of the ray-triangle intersection problem from O(N ) to O(log(N )) at best, N being the number of triangles of the scene. In this thesis context, we have to be even more careful on this aspect since ray queries can represent up to 95% of the computing time on the GPU for a complex scene (containing several millions of triangles). The last 5% is taken by sampling and shading. We review the four main SADS in Section 4.1, to identigy the most tted for our path tracer. Then in Section 4.2, we compare kd-tree and BVH to further argument the choice we made in our path tracer. After that, we describe how we implement our traversal algorithms on the GPU (cf. Section 4.3). Finally, we propose a solution for some noticeable performance gain by encoding a part of the traversal algorithm using a roped-BVH in Section 4.3.2. The uniform grid was introduced by Fujimoto and Iwata [1985]: it subdivides space regularly, having a unique size of cells. This is the simplest SADS. It is easy to implement and the number of cells is fully controlled by the user, by specifying the size of the cell and knowing the 3D scene extent.

Intersection

Intersecting a uniform grid is pretty straightforward, one just needs to intersect each cell along the ray in a "closest to farther" cell order. Once an intersecting cell is found, all triangles belonging to that cell are tested for intersection. The uniform grid accelerates ray tracing performance by avoiding intersection test with all triangles contained in cells not intersected by the ray, or farther than the closest triangle intersected.

Octree

The octree was introduced simultaneously by two researchers, [START_REF] Hunter | Ecient computation and data structures for graphics[END_REF] and [START_REF] Reddy | Representation of three-dimensional objects[END_REF]. It is similar to the uniform grid but adds the potential of having dierent cell sizes, reducing the empty space subdivision seen with the regular grid. The process of creating an octree is similar to the uniform grid, it starts by dening a starting cell size, then it subdivides each nonempty cell until it reaches a minimum cell size threshold or a minimum number of triangles contained in the cell threshold. The octree cells are stored in an octonary tree, each cell containing 0 or 8 child nodes. A subdivided cell is always split at the middle of each of its three axes X Y Z.

Intersection

The intersection process for an octree is the same as that for a uniform grid, adding the step of going down in the hierarchy of cells when the ray intersects a cell that has been subdivided by the octree creation process.

KD-Tree

Kd-trees were introduced by [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF]. To date, they are, with BVH the most used acceleration data structures due to good ray tracing performances. We give a comparison between these two structures, based on the article by [START_REF] Vinkler | Performance comparison of bounding volume hierarchies and kd-trees for gpu ray tracing[END_REF] in Section 4.2. However, they suer from a high construction cost, due to the diculty of nding the split planes that subdivide voxels, and a high memory consumption.

A kd-tree is, like an octree, a space partitioning structure, but stored as a binary tree. Each node of the tree can contain 0 or 2 child nodes. A kd-tree divides the space with planes perpendicular to one of the coordinate system axes X Y Z. A split plane might cut some triangles, so a kd-tree needs to store multiple references to triangles that fall in multiple nodes. This leads to a potential high memory consumption.

Construction

Spatial median split plane. A simple construction method for a kd-tree is to choose alternatively each axis for the split plane, and set the split position at the median of the current node. This naive construction scheme referred as spatial median splitting does not produce the best trees but is often chosen for its low computation cost.

The Surface Area Heuristic. Many construction heuristics have been introduced to build kd-trees: most of them rely on a metric introduced in 1990 by [START_REF] Macdonald | Heuristics for ray tracing using space subdivision[END_REF], the Surface Area Heuristic (SAH). They make the observation that the number of rays likely to intersect a convex object is roughly proportional to its surface area, assuming that ray origins and directions are uniformly distributed in 3D space, which is the case in our ray tracing scenario. The SAH is computed as follows: An optimal kd-tree would be the one that minimizes the SAH for a scene. However it is impossible to compute all the possible trees for a scene, so, most kd-tree construction algorithms rely on a local and greedy approximation that also acts as a termination criteria, computed as follows:

SAH = n∈N odes SA(n) SA(root) × C i + l∈leaf nodes SA(l) SA(root) × C t ( 4
SAH local = C t × T node -C i + C t × SA(Lef t) × T lef t SA(node) + SA(Right) × T right SA(node) (4.
2) where Lef t and Right are the potential left and right child of the node if it is split, and with:

T node : number of triangles in the node T lef t : number of triangles in the left child node T right : number of triangles in the right child node A positive SAH local means that the intersection cost of all triangles in the node is superior to the cost of intersecting a sub-tree composed of the left and right children of this node, and so the node must be split. Most of the SAHbased kd-trees also add a minimum number of triangles contained in a node, so a node must fulll two conditions to be split: having an SAH local > 0 and enough triangles. The C i and C t are empirical values that are implementation dependent.

Even though, the local SAH is not sucient to construct an optimal kd-tree eciently, because at the construction of each node, several split planes that generate dierent potential sub-tree and so dierent local SAH costs exist. Finding the best split plane in an aordable computing time is not trivial and we do not address this problem in this thesis. A good solution has been introduced by [START_REF] Wald | On building fast kd-trees for ray tracing, and on doing that in o(n log n)[END_REF]. They start with a naive approach that builds a kd-tree in O(N 2 ), then rened it to a O(N log 2 N ) complexity by sorting the triangles at each recursive call of the kd-tree construction. Finally, by using a presorting step, they reach O(N log N ) complexity which is the asymptotic lower bound. It is proven easily that O(N log N ) is the lower bound for a SAH kd-tree construction. Let N be the number of triangles and T(N) the cost of building a tree composed of N triangles. If we assume that a split plane divides N triangles into two bins of N 2 triangles, and that the cost of nding such a plane is in O(N), the equation below proves this lower bound.

T (N ) = N + 2T N 2 + ... + 2 log(N ) T N 2 log(N ) = log(N ) i=1 2 i N 2 i = O(N log(N )) (4.
3) However, the fast construction of an ecient kd-tree is not trivial. A solution has been proposed by [START_REF] Wu | Sah kd-tree construction on gpu[END_REF]. Their method can produce high-quality kd-trees using an SAH-based construction, and is optimized to run on the GPU. However, the construction time with their method is still too slow to consider reconstructing the kd-tree at each frame.

Overview of Spatial Acceleration Data Structures

Traversal on GPU

In this thesis and in our path tracing implementation we focus on BVH, see Section 4.1.5. We did not design any particular algorithm for kd-tree traversal on GPU. For a stackless kd-tree traversal on GPU, one can take a look at the work of [START_REF] Popov | Stackless kd-tree traversal for high performance gpu ray tracing[END_REF], in which they adapt an existing stackless traversal algorithm for kd-tree to the GPU and give a CUDA implementation that achieves good performance.

BSP-Tree

The term BSP-Tree refers to Binary Space Partitioning tree. BSP-Trees generalize kd-trees by adding the capability of the split plane to be oriented arbitrarily. Due to their complex topology and high combinatorial construction process, they are not really practical for ray tracing applications on the GPU.

BVH

All SADS previously introduced in this document are space partitioning structures. Another class of SADS is object partitioning structures: their construction process relies on a hierarchical object representation. Bounding Volume Hierarchy (BVH) belongs to this family.

A BVH is stored as an N-ary tree, each node containing 0 or N pointers to its children nodes. In this thesis, we focus on binary tree BVHs. Each node that does not have any child is a leaf node. To minimize the impact on memory, we store references to triangles only in the leaf nodes.

Memory and node boundaries

Since BVHs do not split triangles, except for some peculiar, more advanced BVH that we will introduce later, they have the property of being bounded in the number of nodes and depth. Furthermore, for a BVH containing N triangles, the worst case is a "left or right comb": a BVH having all its leaf on left or right, as shown in Figure 4.5. For a binary BVH containing N triangles, its maximum depth is N and its maximum number of nodes is 2N -1. This nice property is more than appreciable for a GPU path tracer, where, currently, memory is more limited than on a CPU.

Construction

To construct our BVH, we rely on the SAH introduced in Section 4.1.3. We x the minimum of triangles in a leaf to six, C t = 1.0 and C i = 2.4. As these parameters really depend on the implementation of the BVH traversal, a lot of dierent values can be found in the literature. For instance [START_REF] Aila | Understanding the eciency of ray traversal on gpus[END_REF]] use a minimum of eight primitives with C i = 1.2 and C t = 1.0, whereas [blender.org, 2017] uses a minimum of one primitive with C i = 1.0 and C t = 1.0.

We do not build a unique BVH for the whole 3D scene, but a BVH for each dierent geometry present in the scene. We do take care of multiple instances of geometries, allowing them to reference a unique BVH. By doing so we can easily move objects in the scene without any reconstruction of the SADS. This method has some drawbacks, as it leaves us with as many BVHs as dierent geometries in the 3D scene. To counter that, we could build a unique BVH over the whole set of triangles of the 3Dscene. It is equivalent as considering the 3D scene as a unique geometry. But, this would force us to add a reference to a material for each triangle. By having dierent geometries, and using an instance system, we can assign a material to a whole object. An object is a set of a material, a model matrix and a reference to a geometry.

Meta-BVH

We tried to add a meta-BVH over this at hierarchy of BVHs, to better organize the scene. Our meta-BVH is built using the same heuristic than our BVH. Instead of containing triangles its nodes reference object BVHs. Unfortunately, as presented in Table 4.1, results were not conclusive. The advantage of our meta-BVH was noticeable only in pathological cases like the multiple dragon scene presented in This is mostly due to the repartition of objects in the scene. In the San Miguel scene for instance, there are 253 objects which have bounding boxes with an extent that covers the whole scene, leading to 253 huge bounding boxes. Obviously, adding a Meta-BVH to such a chaotic scene hierarchy does not improve performance.

In conclusion, this Meta-BVH heuristic, to be ecient, needs an algorithm that reorganizes the scene hierarchy and creates or merges objects.

Faster Intersection: Roped-BVH

To improve Algorithm 4.3, we have made the choice of encoding the traversal in the tree. In fact looking at Algorithm 4.3 and Figure 4.7 we see that the depth traversal of the tree is deterministic and that case B in Figure 4.7, which illustrates the goToNextNode() function, can be very expensive. Thus a simple solution is to store pointers to next node to test and replace the goToNextNode() function by only one node fetch (see Algorithm 4.4). This new pointer hierarchy is illustrated in Figure 4.8, where the orange curves represent the new pointers. This new hierarchy has another big improvement, it consumes less memory, as there is no need to store a pointer to the parent anymore. The left and right child pointers are replaced by two pointers: one to the next node to test if the ray intersects the node, the other to the next node to test if the ray does not intersect the node. In the case of a leaf node, only one pointer is needed, the next node to test. This new traversal algorithm gives a better performance on GPU. We did see a 1.3 to 1.9 acceleration factor (depending on the scene), compared to the traversing Algorithm 4.3. A benchmark is presented in Table 4.2. The same approach has been recently published by Torres et al. [2009], without receiving much attention.

Roped-BVH Memory Layout on the GPU

We present here in details the memory layout of our roped-BVH. It is stored in two GPU buers: the node buer and the BBOX buer. The node buer stores the tree structure (i.e., the links between each node of the tree and the node description). This node buer is a 1D array of integers. A link in the tree is just an integer indicating in which cell the pointed node is stored: this is equivalent to a pointer mechanism. A non-leaf node is then stored on four integers:

• a pointer to the next node to traverse in case the ray intersects the current node

• a pointer to the next node to traverse if the ray does not intersect the current node

• a pointer to the BBOX of the node in the BBOX buer

• the number of triangles in the sub-tree starting at this node.

The number of triangles is stored as a negative value do dierentiate non-leaf node and leaf node. The two pointers to the next nodes correspond to the orange links in Figure 4.8. In the case of a leaf node, we add to this four integers a set of pointers to the triangles contained in that node. Finally, the BBOX buer stores as oating point values the XYZ coordinates of the bottom-left and up-right corners of each BBOX of each node of the tree.

As future work, it would be interesting to see if we can gain even more performance by replacing the integer (4 bytes) pointers by smaller data types like short integers (2 bytes). This, of course, can be done only if the number of nodes in the tree does not exceed the maximum value that can be represented by this smaller data type. By doing so it potentially would reduce both memory consumption and traversal time because a node fetch would probably be faster.

Conclusion and Research Perspectives

We have presented in this chapter the four main SADS, and described how we use the BVH that we have chosen for our path tracer. Even though SADS were not the main area of study in this thesis we tried to get the best performance from the BVH, for instance by developing the roped BVH. To push further we could have used a more advanced SADS, like the SVBH from Stich et al. [2009]. They describe a hybrid SADS, that adds to the BVH the possibility to split triangles. This SBVH is somehow similar to a kd-tree and is probably, up to now, one of the most powerful acceleration data structures. Another potential good solution for a robust production renderer might be to switch between a kd-tree and a BVH depending on the 3D scene, and 80 Arthur Dufay this is the solution NVIDIA took in Optix [START_REF] Parker | Optix: A general purpose ray tracing engine[END_REF]). In that sense, some research on how to choose a more appropriate SADS might be more than helpful. In addition to that, maybe the SAH heuristic needs to be revisited: a new heuristic, even empirically developed from benchmarks of ray tracing on 3D scenes on the newest GPU hardware might be giving good results in constructing SADS. Indeed, as previously explained in Chapter 3, dierent GPU architectures induce dierent memory layouts and performances: one optimization that works on a GPU might not be as good on another, so an adaptive solution might be the best option.

Once again one of the main limitations of the GPU is memory cache. We will see in Chapter 5 how we further enhance performance by increasing cache coherency.

Chapter 5

Random Number Generation on the GPU

In the rst chapter, we saw that computing realistic images can be achieved by using Monte-Carlo integration. However, this requires generating random numbers that will be used as samples for the Monte-Carlo integrator. Generating random numbers on the GPU is not trivial, and we have to take good care of the properties of the generated samples as well as the computation time involved by the random number generator (RNG). We also pointed out in the previous chapters that cache coherency is crucial to obtain maximum performance. Random samples are used to produce random ray directions and thus sample generation and ray tracing performance are tightly linked.

In this chapter, we present our work to obtain more coherence between rays. We rst introduce a mathematical tool to evaluate the quality of a sampling (cf. Section 5.1). Then in Section 5.2, we describe several methods to generate random numbers on the GPU and see their benets when they are used for a GPU path tracer. After that in Section 5.4, we present our new decorrelation method. It depends on some parameters, relative to the sampling used, and can be adjusted semi-automatically.

Discrepancy

Discrepancy

In Chapter 1 we saw that importance sampling can be used to improve eciency of a Monte-Carlo integrator. However, when no particular information is known about the function of which we try to estimate the integrand, a uniform sampling is done. For instance, this is the case when we try to solve the rendering equation for a perfectly diuse surface. The lambertian surface reects light equally in all directions, as we have seen in Figure 2.6 and no information of the visibility is known prior to rendering. In that case, the more the samples are well distributed the more ecient the Monte-Carlo integrator will be on average.

To ensure that samples are well distributed, a measure known as the discrepancy, and noted D N (P ), has been introduced. It is well presented by [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]. Intuitively, the discrepancy is the biggest dierence between the density of samples expected and measured in a sub-volume J; it is computed as follows:

D N (P ) = sup B ∈ J A(B; P ) N -λ s (B) (5.1)
where P is the sample set, N = card(P ), A(B; P ) are all the samples of P that fall in the sub-volume B and λ s (B) the S-dimensional Lebesgue measure. J is the set of S-dimensional intervals of the form:

J : s i=1 [a i , b i ) = {x ∈ R s : a i x i < b i } with 0 a i < b i 1 (5.2)
Another measure, known as the star discrepancy, is more practical. In fact it is computed on a more restrictive set of intervals J * :

J * : s i=1 [0, b i ) = {x ∈ R s : 0 < b i } with 0 < u i 1 (5.3)
Star discrepancy corresponds to consider only sub-volumes of the integration domain that start with a vertex at the origin [0] s . As the star discrepancy is much easier to compute we will use it rather than the discrepancy to assess the quality of the sampling used. Since the star discrepancy of a sampling guarantees that it is well distributed, when estimating the integrand of a function with it, it ensures that we will not miss important values of the function.

5.2

Random Number Generation on the GPU

Fast Pixel-based Techniques

On the CPU there are built-in functions such as C-language rand() from the STL or libC, or the more recent C++11 mt19937, an implementation of the 84 Arthur Dufay

Random Number Generation on the GPU

Mersenne twister algorithm from [START_REF] Matsumoto | Mersenne twister: A 623-dimensionally equidistributed uniform pseudo-random number generator[END_REF]. On the contrary, no random functions are provided on the GPU. Thus to generate random numbers we often rely on a combination of trigonometric functions and permutation of bits. These techniques could be classied in pixel-based techniques because they rely on a seed, often determined with the pixel-ID to generate a random number. For instance, one could use the method described in Algorithm 5.1. This technique has the advantage of being easy to implement, predictive and has a low computational cost.

Algorithm 5.1: A pixel based RNG, with its associated seed function.

Wang Hash

Another interesting RNG was introduced by Wang [1997] and is computed as presented in Algorithm 5.2. Algorithm 5.2: A pixel-based RNG, with its associated seed function.

Low Discrepancy Sequences

Low discrepancy sequences (LDS) are an alternative to fully random numbers for Monte-Carlo integration. Their intrinsic property is that for any sub-sequence S : x 1 , ..., x N , S has a low discrepancy. Thus, using an LDS to compute a Monte-Carlo integrator improves eciency as explained in Section 5.1. A method that uses an LDS instead of a fully random sequence to 86 Arthur Dufay 5. Random Number Generation on the GPU compute a Monte-Carlo integrator is called a quasi-Monte-Carlo method. We expose in this section three well-known low discrepancy sequences.

Van der Corput Sequence

We rst introduce the van der [START_REF] Van Der Corput | Verteilungsfunktionen. I. Mitt[END_REF] sequence, that is dened as:

g b (n) = L-1 k=0 d k (n)b -k-1 (5.4)
where b is the base in which number n is represented.

Hammersley Sequence

The [START_REF] Hammersley | Monte carlo methods for solving multivariable problems[END_REF] sequence is an LDS based on the van der Corput sequence. It is computed as follows:

x(n) = g b 1 (n), ..., g b s-1 (n), n N (5.5)
where b 1 , ..., b s-1 are co-prime integers greater than 1 and N is the number of samples. We see here that using a uniform distribution as the rst dimension for the Hammersley distribution restricts its usage. Indeed the number of samples has to be decided prior to rendering and all the samples must be taken. Otherwise, the Monte-Carlo integrator will be biased if the pdf is not corrected accordingly to the non-taken samples. Due to this restriction, we did not use this sequence in our path tracing. However, the Hammersley sequence is still a good choice when the number of samples is xed: it is well distributed and really easy to compute.

Halton Sequence

The [START_REF] Halton | On the eciency of certain quasi-random sequences of points in evaluating multi-dimensional integrals[END_REF] sequence solves the problem of the xed number of points in the Hammersley sequence by using only numbers from the van der Corput sequence. It is computed as follows:

x(n) = (g b 1 (n), ..., g bs (n)) (5.6)
where b 1 , ..., b s-1 are co-prime integers greater than 1. It is well distributed and easy to generate on a GPU, but as shown in Figure 5.3, the sampling quality decreases when going into higher dimensions.

Sobol Sequence

The Sobol sequence is another LDS introduced by [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. Before describing its properties we must introduce two denitions from [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF], the (t,m,s)-nets and the (t,s)-sequences. In other words, in our case, a (0,m,2)-net will hold a unique point in each sub-interval of side length b -m of the integration domain. This property guarantees a low discrepancy value, and so a fast convergence in a GPU path tracer.

The Sobol(0,2) sequence is a (0,2)-sequence in base 2, it implies that each successive set of 2 m points is a (0,m,2)-net, giving a fast convergence when using it in a path tracing engine.

However, due to its complexity to generate it, it is tedious to use it in a GPU path tracer. Its values must be precomputed on the CPU and uploaded to the GPU. It also has the same problem than the Halton sequence: its distribution quality decreases as we get to higher dimension (cf. Figure 5.4).

Decorrelation

Introduction

In the context of path tracing, when using low discrepancy sequences we have to deal with decorrelation of samples between the dierent pixels of the generated 88 Arthur Dufay This is why we did not consider it as a solution in our context of previsualization. However, some work has been done by [START_REF] Sadeghi | Coherent path tracing[END_REF] to minimize those coherency artifacts. By interleaving several coherent sequences they break the structured noise. To sum up, to use an LDS in a path tracing engine we have to decorrelate samples.

Decorrelation Techniques

We now present some decorrelation techniques useful in the context of path tracing. Rotation matrix M can be computed using one of the pixel-based techniques described in Section 5.2.1. The main advantage of this method is that it preserves the discrepancy property on the disc. On the other hand, this decorrelation technique does not cover properly the integration domain. As shown in Figure 5.7, we observe undesirable circular patterns and areas that never get sampled.

90

Arthur Dufay Let M the number of pixels in the image, N the number of samples in S. The decorrelated sub-sequence S p i for pixel p i is then written as follows:

O : [u 1 , ..., u M ], u i ∈ N S : [x 1 , ..., x N ] S p 1 : [x u 1 +1 , ..., x u 1 +N ]...S p N : [x u N +1 , ..., x u M +N ] (5.7)
This technique still has some drawbacks. When using it on a GPU path tracer, several threads working in parallel in the same warp ask for samples that might fall far from each others in the LDS. In the case of the Halton sequence, for instance, generating such samples with dierent IDs in the LDS in parallel leads to code divergence. In other words, this decorrelation technique can potentially increase the computing time of the RNG. In our opinion, this technique is also not sucient to fully get rid of coherency artifacts. As many pixels potentially use the same or a close oset, they will have too many samples in common, which introduces some bias in the Monte-Carlo estimator.

Cranley Patterson Rotation

Cranley Patterson rotation, introduced by [START_REF] Cranley | Randomization of number theoretic methods for multiple integration[END_REF], is a well-known decorrelation technique. It uses a D-dimension vector δ of ran-dom values to alterate a sampling of dimension D. This random vector is chosen once and for all the samples. Applying Cranley Patterson rotation to a sampling S to obtain a new sampling S is written as follows:

δ : [u 1 , ..., u D ], u i ∈ R S : [x 1 , ..., x N ] S : [x 1 + δ -x 1 + δ , ..., x N + δ -x N + δ ] (5.8)
The method is called a rotation because a oating modulo is applied on the altered value, as shown in Figure 5.8. We found it to be the most useful technique to generate quickly new samples on the GPU. One just needs to precompute on the CPU a buer of random values with as many values as pixels in the image to render. Then, when using the same RNG for each pixel, for instance, a Halton sequence, one just needs to alterate the output of the RNG with the value stored in the Cranley Patterson buer.

5.4

A GPU Cache Friendly Decorrelation Technique -Micro Jitter

Motivation

We saw in the previous section that Cranley Patterson rotation can be used to decorrelate samples in path tracing. However, using it straight as it is on the GPU could have some impact on the performance. Actually, if we compare it to coherent path tracing, as presented in Section 5.3.1, it tends to generate rays that are fully incoherent (cf. Figure 5.9).

Arthur Dufay

Figure 5.9: The Cranley Patterson rotation method induces totally random numbers and so fully incoherent rays. This breaks the GPU cache coherency as paths are computed in parallel.

To solve this problem, one can rely on coherent path tracing (i.e., avoid any decorrelation method). Indeed it greatly helps in preserving the GPU cache coherency as shown in Figure 5.10. However, as presented in Section 5.3.1, it is not a solution for previsualization because, until the image is fully converged, it exhibits some structured noise. Another solution is to rely on ray-reordering techniques, but their implementation on the GPU is quite tedious.

Seeing the good performance that coherent path tracing can give, we wanted to nd a solution in between coherent path tracing and Cranley Patterson rotation, i.e., a solution that can remove the structured noise of coherent path tracing, but at the same time, keep some GPU cache coherency. We came up with the idea of a micro jitter method that we describe in the next section.

Method Description

Our micro jitter method works similarly as Cranley Patterson rotation. It is a decorrelation technique. The only dierence with Cranley Patterson is that the jitter vector δ, used to alterate the sampling, is randomly chosen in a well-dened range. By controlling the amount of jitter we apply, the new set of samples generated with our method tends to generate rays that are similar in 3D space. This similarity between rays greatly helps in reducing cache misses when rays are traced in parallel on the GPU (see Figure 5.12). Quasi-similar rays that will be traced by threads belonging to the same warp will likely hit the same object in 3D space (i.e., they will fetch the same BVH nodes) and will be faster to trace than totally incoherent rays. This is especially true for rst bounces of paths.

Adjusting the Jitter Radius

General Formulation To apply our method, the amount of jitter must be selected carefully. As shown in Figure 5.13, if the amount of jitter is too small, by superposing several jittered sets of samples, we clearly see that some parts of the integration domain are not covered. Indeed, a jitter radius too small generates new points in a too small region. If the average space between samples is larger than this region, our method cannot generate points that fully cover the integration domain. The jitter radius must be chosen depending on the distribution of the input sample set.

However, the amount of jitter must be as small as possible to maintain a protable gain in performance. As the amount of jitter gets bigger our 94 Arthur Dufay Figure 5.12: Our micro jitter technique preserves cache coherency by generating similar rays in 3D space, while avoiding artifacts of coherent path tracing.

Figure 5.13: Our micro jitter applied once (left), 150 times with a too small jitter (middle), with a proper amount (right). The original distribution is shown in blue.

A GPU Cache Friendly Decorrelation Technique -Micro Jitter

method tends to reproduce the Cranley Patterson rotation behavior, and so, the ray tracing performance drops. On the other hand, a jitter too small would reproduce the structured noise of coherent path tracer.

A perfectly distributed sampling will leave, on average, an empty hypercube around each sample of size N -1 S , N being the number of samples, s the number of dimensions. For a 2D sampling it would be 1 √ N . Our jitter radius µ must then be set to: µ = K √ N with K = 1 for a perfectly distributed sampling. This can be generalized to any dimension with the formula: µ = KN -1 S . K depends on the quality of the sample set.

Application to the Halton Sequence However, when dealing with LDS as the Halton sequence, the distribution is not perfect. In consequence, we must set µ to a larger value. As the star discrepancy is a measure of the distribution of samples, we use it to select our amount of jitter. We nd out empirically that the star discrepancy of the Halton sequence is roughly proportional to f (x) = 2.5 √ N (cf. Figure 5.14). We could have set our jitter radius to µ = 2.5 √ N but, in practice, we do not have any information about the shape of the empty space between samples. They can be perfectly squared or very elongated. Therefore, to ensure that we do not miss part of the integration domain, we double this jitter radius. For the Halton sequence the samples are then jittered by a vector: δ ∈ [-µ, µ] with µ = 2.5 √ N . To use our micro jitter method with another sample set, one just needs to adjust the jitter radius with the star discrepancy of the desired sample set.

Results

Performance Analysis

As shown in Figure 5.17, we tested our micro jitter method over several LDS and obtained a better performance. The results average the experiment conducted on three dierent scenes (all shown in Figure 5.16): the interior scene (85K triangles), the museum scene (15M triangles) and the dragon scene (870K triangles). We tested our technique with ambient occlusion (AO) computation and path tracing with a maximum of three indirect bounces. A close-up view of path tracing using our method is shown in Figure 5.18.

As previously explained, the amount of jitter is linked directly to the number of samples in the sequence. As the sampling count increases, the amount of jitter gets reduced and performance increases. The upper bound for this performance gain is the performance of coherent path tracing.

We can also notice that the performance gain is better for AO computation than path tracing (PT) computation. Indeed, AO computation requires launching only one secondary ray, not a full path, whereas PT launches one to three bounces in our experiment. As our method increases performance by taking prot of ray coherence, it is quite straightforward that the rst bounce of the path gets a better performance gain, hence the better performance gain for AO computation. If we carried out tests with more than three bounces for PT, we would have seen a smaller performance gain with our micro jitter method.

Figure 5.17: Performance gains of our method on several rendering algorithms: Ambient occlusion (AO) and path tracing, tested on dierent samplings.

Error Analysis

We also conducted some analysis on the images generated with our method. To do so, we computed the RMSE (Root Mean Square Error) of image dierences of ambient occlusion images computed with our decorrelation method and a reference image. We did the same thing with images computed with the

Limitations

We presented here a method that improves cache coherency in path tracing. However, our method has some drawbacks. Firstly, with our micro jitter, the number of samples must be determined in advance to adapt the jitter radius (cf. Section 5.4.2). This is a huge limitation, but it can be overcome, using batches of samples. Secondly, our method is scene topology dependent. Incoherent scenes with small triangles, like the Hairball scene (cf. Figure 5.20), have too many variations of triangle normals from one pixel to another. This generates totally random directions for the secondary rays, whatever the decorrelation technique chosen.

Figure 5.20: The Hairball scene (2,880,000 triangles) does not benet from our method. Too many dierent triangle normals generate totally random rays in any cases.

Conclusion and Future Work

We presented here a simple way to improve ray tracing performance in a GPU path tracer. The development of our method led to the submission of a patent (Micro-Jittering for GPU-Friendly Monte Carlo Multi-Dimensional Integration Problems) and a publication in a major conference [START_REF] Dufay | Cache-friendly micro-jittered sampling[END_REF]. It was successfully tested in the open source path tracing engine Blender Cycles (see blender.org [2017]), with just a few lines of code. Our method is indeed much easier to implement than a full ray sorting algorithm like [START_REF] Garanzha | Fast ray sorting and breadth-rst packet traversal for gpu ray tracing[END_REF]. This contribution has been well received by the scientic community and has been added to the Blender Cycles roadmap for the next release. As future work, we would like to compare the benet of our method to a full ray sorting solution. Another point would be to address more in detail how to adjust automatically the jitter size to any sample distribution, especially in the case of adaptive sampling when the number of samples is not xed.

102

Arthur Dufay

Conclusion

In this thesis we addressed the previsualization of VFX. Our main goal was then to enhance an already existing 3D rendering platform with global illumination capability. For that we chose to focus our interest on the path tracing algorithm. Its implementation on dedicated graphic hardware can, indeed, truly help VFX artists in their designing tasks, by giving them a faster feedback and a resulting image closer to the nal render algorithm output. However, the setup of such rendering engine in an industrial context was not trivial. For that, we had to tackle several aspects of the rendering engine: BRDFs and materials (cf. Section 1.2.3), GPGPU (cf. Section 2.2.1 and Chapter 3), spatial acceleration data structures (cf. Chapter 4) and random number generation (cf. Chapter 5). This lead to several contributions that are summarized here.

Contributions Summary

First of all a full path tracing engine was successfully implemented inside the 3DCast platform. Inside that, a solution to increase interactivity was proposed, our quad-tree computation (cf. Section 2.2.4).

A rst patent on tracing shadow rays in ray tracing applications was submitted (cf. Section 3.6). It clusteres shadow rays by light sources and launches rays in a reverse manner (from light sources to surfaces) to improve ray tracing performance.

An optimization of the traversal algorithm of a BVH on GPU was proposed in Section 4.3.2. Using an encoded traversal inside the spatial acceleration data structure it further increases the eciency of ray tracing intersection test in our rendering engine.

This thesis also lead to the development of a new decorrelation technique that proves to be protable for several multi-dimensional integration problems including path tracing and SSAO computation. This contribution has been published in a major conference [START_REF] Dufay | Cache-friendly micro-jittered sampling[END_REF], and submitted as a patent: Micro-Jittering for GPU Friendly Monte Carlo.

Finally, a third patent, to increase speed of shadow computation is presented in Annex B: it is submitted.

Implementations

Several topics of our path tracing engine were not completed during this thesis. This leads to several future projects.

First, the support of more complex material denitions, such as BSSRDF (e.g., skin behavior), is essential to deal with complex surface description, as found in a production renderer.

In pair with such materials we would like to implement another path tracing GPU pipeline featuring the wavefront idea of Laine et al. [2013]. That would have been benecial only when dealing with complex materials. For spatial acceleration data structure, we use a quite standard BVH. The implementation of a more advanced one, using some splitting heuristic, that increases BVH performance and makes it closer to a kd-tree would be an interesting point to investigate.

Finally, even though the exibility of an easy setup oered by OpenGL Compute Shaders was protable, we would like to implement our solution using the NVIDIA CUDA framework. We believe that the proling and debugging tools accessible for this GPGPU framework would be a tremendous asset in our case.

Research Topics

Bidirectional solution The major research area we would focus on next is bidirectional algorithms such as bidirectional path tracing (BDPT), or even a more complex solution, that combines both path tracing (PT) and BDPT. We thought of using heuristics to decide whether or not a camera path must try to connect to a light path. One of these heuristics could be based the material encounter along the path. Remember that a specular path might be dicult to connect or even impossible to connect with highly specular BRDFs (cf. Figure 1 .11). Another heuristic might be the quantity of energy transported along the camera path (i.e., only low energy camera paths might need connections to light paths). Using such algorithm we could potentially greatly improve performance of a BDPT by making some paths bidirectional and others unidirectional. Indeed, coupling two buers to sort bidirectional paths and unidirectional paths, we could save some computation time.

Hybrid algorithms Another point we would like to investigate is hybrid algorithms. We think that the many lights methods truly have a potential for VFX previsualization. Combining these techniques with path tracing to remove their artifacts (cf. We started some test at the end of this PhD where we used some heuristics to decide for each pixel of the image whether to use VPLs or path tracing for indirect illumination. We used a heuristic based on the harmonic distance. We precomputed it using some ray tracing method. Then indirect illumination is computed using PT if the harmonic distance is under a desired threshold, otherwise VPLs are used (cf. Figure 5.22).

The next step would be to nd a way to have a smooth transition between the two methods. that are lit or shadowed. For remaining unclassied view samples they further perform an intersection test against the rasterized triangle by retrieving its coordinates through its ID. For remaining pixels classied as uncertain they cast shadow rays using a GPU ray tracer to solve their lit status.

Drawbacks Their method requires triangle indexing of the scene and storage of the triangles in a geometry buer. If vertices are animated, this buer needs to be updated and may introduce some computational overhead. Moreover many dierent triangle IDs can be stored locally in a shadow map especially in presence of unorganized geometry (trees, leaves for instance). Storing unorganized triangle IDs lead to scattered memory accesses that are proved to be inecient on GPU architectures. We suspect this as the bottleneck of their method because the reported speed up does not show so much improvement in rendering times. Best speed up reported is 1.46× over a full ray tracing solution. The use of the triangle coverage ag is prone to errors for slanted surfaces.

Finally their technique performs excessive ray tracing computations in areas lying along shadow edges. Computational overhead becomes prohibitive in presence of models with geometry density similar or slightly greater than the shadow map resolution. We propose a simple and ecient solution to handle this (cf.

B.2.2 Selective Ray Tracing

Lauterbach et al. [2009] use also a conservative shadow map to accelerate shadow ray tracing computations. They propose to analyze the 8 depth samples around a depth sample to detect the presence of shadow edges. It is done by comparing the maximum absolute depth dierence with the minimum depth variation determined by the far and near planes. Depending on the surface inclination regarding the light source they determine pixels that need to be ray traced or not.

Drawbacks The method fails in detecting small holes in geometry due to the nite resolution of the shadow map and the lack of geometric information. Moreover it is sensitive to shadow map bias.

B.3

Technical Problem Solved by the Invention

The proposed invention takes advantage of ecient GPU rasterization engines combined with GPU ray tracing techniques to drastically accelerate rendering times of ray traced shadows. It relies on a fast and precise lit/shadowed pixel 112 Arthur Dufay classication using conservative shadow maps, that outperforms previous solutions. Ray tracing shadow computations are then limited to a small subset of pixels that signicantly reduces ray tracing operations over previous solutions. Our solution does not introduce any artifact and could achieve real-time performance in complex scenes at quality equivalent to full ray traced shadows.

B.4 Proposed Solution

Our solution consists in rasterizing the shadow map in a conservative way and explicitly storing the entire triangle denition in shadow map pixels. By proceeding like this we avoid scattered random accesses into GPU memory.

Secondly we propose a modication of depth fragment generation that improves the lit classication and gives a better estimate of maximum ray length for ray traced shadows. Finally by considering neighborhood of shadow map pixels we further improve the classication lit/shadowed of pixels at shadow edges. GPU ray tracing operations are then limited to a small subset of pixels in the image resulting in drastically reduced computation times.

B.4.1 Overview

Our solution works at three stages:

• In the rst stage we render the conservative shadow map with explicit triangle storage using a modied minimal conservative depth.

• In the second stage, the scene is rendered. Using a fragment shader, we proceed to pixel classication by querying the shadow map.

• In the third stage, GPU ray tracing of shadow rays is performed on unclassied pixels with limited ray traversal distance.

B.4.2 Conservative Shadow Maps with Explicit Triangle Storage

We describe in this section the generation of the conservative shadow map. We briey recall the principles of conservative rasterization, then we describe how we pack the entire triangle denition in a single RGBA shadow map pixel.

Then we describe a new method that computes a minimum depth value at a pixel that is more consistent regarding the triangle coverage within the pixel.

Conservative Rasterization

The rst pass consists in generating the conservative shadow map. Standard rasterization evaluates a triangle at pixel centers. A triangle that intersects a B.4. Proposed Solution pixel but does not overlap its center would not produce any fragment. Consequently no information will be written in that pixel. Conservative rasterization guarantees that fragments will be produced for any closer triangle in the pixel area. Techniques to perform conservative rasterization have been proposed by [START_REF] Hasselgren | Conservative Rasterization[END_REF] and Hertel et al. [2009]. They both rely on the same idea: move slightly triangle edges forward in their normal direction by a length of half a pixel width.

For shadows, conservative rasterization of shadow maps allows rapid classication of lit area. It guarantees that depth of the closest visible triangle will be written whatever the triangle coverage inside a pixel. Therefore if the depth of the view sample projected in light space is less or equal to the depth stored in the pixel, it guarantees that the pixel is fully lit.

Explicit Triangle Storage

At the dierence of Hertel et al. The idea consists in storing for each pixel of an RGBA shadow map the 2D coordinates of the projected triangle, its depth evaluated at the pixel center and the compressed partial derivatives in shadow map space. Using this information, we are able to perform a shadow ray intersection directly against the triangle stored in the pixel. The main advantage is that it requires a single texture fetch to access the triangle information instead of deferred random accesses using triangle IDs.

Minimum Depth Fragment Generation

When performing conservative rasterization we have to choose carefully which depth value to write in the depth buer. [START_REF] Hertel | A hybrid gpu rendering pipeline for alias-free hard shadows[END_REF] use a slope scale depth correction that shifts triangles towards the light source such as the depth value at pixel center corresponds to the minimal value found in the pixel area crossed by the triangle plane. The choice of minimum depth is used for the classication of lit pixels. Thanks to conservative rasterization, if the depth of the view sample projected in light space is less or equal than the minimum depth stored in the corresponding pixel we conclude that the pixel is fully lit.

However if we consider a surface in front of slanted surfaces, as depicted in Figure B.3, the slope scale depth correction may spawn incorrect depth values. In this situation pixels on the green surface cannot be classied as lit resulting in unnecessary expensive ray tracing operations. To that end we exploit the partial depth derivatives computed at the shadow map generation to speed up nding of the minimum depth. Depth derivatives indicate decreasing depth variation along the shadow map axis. They form also a gradient vector in 2D that indicates the decreasing direction in the pixel. We simply use this vector to select the triangle edge that holds the smallest depth. It is identied as the one with normal vector closest to the gradient vector.

Once identied, we perform the intersection of this edge with the box surrounding the pixel and determine depth at intersection points. Depth at intersection point m is computed as follow:

depth i = p Z + (m x -c x ) d p Z dx + (m y -c y ) d p Z dy (B.2)
• If two intersection points are found, we take the minimum depth of both (cf. • If one intersection is found then we conclude that a triangle vertex is inside the pixel. We then take the minimum depth from the vertex location and the intersection point (cf. • If no intersection is found, the edge is either fully inside or fully outside the area. If the edge is fully inside we take the minimum depth from the two vertices location.

• If edge is fully outside, we then select the adjacent edges and perform intersections with surrounding pixel area box.

• If only one intersection is found on each box side pointed by the negative gradient direction in axis X and axis Y, we then conclude that the corresponding box corner is covered by the triangle and take minimum depth at this corner (cf. • Otherwise, if two intersections are found we take the minimum depth at intersection points and vertex. If four intersections are found, we take the minimum depth among the four intersections (cf. • If no intersection is found, the triangle fully covers the texel area. In this case, the minimum depth is found at texel border using Equation B.3.

Implementation

The generation of our conservative shadow map is implemented using a geometry shader for the computation of triangle conservative expansion and a fragment shader for the triangle encoding and minimum depth computation. To overcome these problems we propose to query the eight pixels in the neighborhood and retrieve the surrounding triangles recorded in the shadow map. Because these triangles may potentially cover the pixel at projection location we perform direct shadow ray intersection on this triangle set. If an intersection is found, the pixel is classied as shadowed avoiding the need of a full ray traced operation.

B.4.4 Full Classication

According to the previous sections, the complete pixel classication works as follows:

1. Query pixel of shadow map at view sample p projection.

2. Compute the minimum depth or retrieve it directly from the depth buer. 4. If not, we test the intersection with the triangle stored in the texels as described in Section B.4.2.

5. If an intersection is found, the view sample is classied as shadowed.

6. If not, we test the intersection with the triangles stored in the eight surrounding pixels.

7. If an intersection is found, the pixel is classied as shadowed.

8. Otherwise the pixel is classied as uncertain.

The full classication can be implemented in a fragment shader at rendering time or in a computing kernel (OpenCL, CUDA or Compute Shader) in image space using a G-Buer.

B.4.5 GPU Ray Tracing of Shadow Rays

For areas where the lit status of view samples remains uncertain, we spawn a shadow ray toward the light source and look for an intersection. If an intersection is found, the pixel is considered in shadow otherwise it is considered as lit.

The shadow ray intersection relies on a GPU ray tracer and spatial acceleration structures for faster intersection determination. Basically a ray traverses a tree of sub-spaces (e.g., octree, kd-tree, LBVH, . . . ) until reaching a leaf that contains geometry. Intersection test is done on each triangle contained in the leaf. [START_REF] Hertel | A hybrid gpu rendering pipeline for alias-free hard shadows[END_REF] propose to limit the length of the ray to prevent visiting empty sub-spaces during the ray traversal and thus improve eciency of the GPU ray tracer. This length is determined by the smallest value stored in • Renders high quality shadows at real-time or interactive framerates.

• Does not produce any artifact compared to previous solutions.

• Performs faster than previous solutions.

• Compatible with any ray tracing or deferred rendering pipelines.

• Easy to implement on graphics hardware.

• Easy to implement on production renderers.

Figure 1 :

 1 Figure 1: Left: CG advertisement: Canal+ The Bear, courtesy of Mikros Image. Right: CG used to predict architecture design, courtesy of NVIDIA.

Figure 3 :

 3 Figure 3: Example of the progress made in the context of computer-animated movies, from Toy Story 1995 (Left) to Finding Dory 2016 (Right). Both courtesy of Pixar Inc.

Figure 1

 1 Figure 1.1: A rendering of a triangle mesh: "The Stanford Bunny" (69K triangles) from Stanford University [1993]. Left: solid view, Center: wire-frame view, Right: closeup of wire-frame view.

  Figure 1.2: A nal render image from Guardians of the Galaxy, courtesy of MPC.

  Figures 1.4 and 1.5) such as global illumination and caustics. A global illumination algorithm generates images that sum the contributions of direct lighting and indirect lighting.

Figure 1 . 3 :

 13 Figure 1.3: Direct lighting shown on the hemisphere centered on A, lit directly by the area light. Indirect lighting shown on the hemisphere centered on B, lit by photons bouncing on the scene.

•Figure 1 . 4 :

 14 Figure 1.4: Path tracing image demonstrating global illumination. Top: direct illumination only. Bottom: direct and indirect illumination.

Figure 1

 1 Figure 1.5: A sketch of a caustic path LSDE (top), and its rendered version, a caustic from a gold ring lit by a point light source showing indirect illumination only (bottom).

Figure 1

 1 Figure 1.6: Left: one steradian, it subtends an area of r 2 , Ω = A r 2 = r 2 r 2 = 1 sr. Right: Representation of the domain of integration Ω for the rendering equation. The Radiant energy noted Q is the energy of electromagnetic radiation. Its unit is in J . The Radiant ux noted Φ is the radiant energy emitted, reected, transmitted or received, per unit time. Its unit is in W or J • s -1 .

  Figure 1.7: Estimating π with Monte-Carlo, 82 points inside the disc, 100 points in total, 82 100 * 4 = 3.24 ≈ π.

1. 4 .

 4 Figure 1.8: Ray tracing principle with a pinhole camera.

  Figure 1.9: Path Tracing (Left) and Light Tracing (Right) act in reverse order.

Figure 1 .

 1 Figure 1.10: Bidirectional Path Tracing constructs paths by connecting eye path and light path.

Figure 1 .

 1 Figure 1.11: Caustic connection problem. Left: a possible caustic connection. Right: a zero energy path. BRDF distributions are shown in orange lobes and hemispheres.

Figure 1 .

 1 Figure 1.14: Back-face culling applied on a cube. Green polygons are kept, red ones are seen from behind and are discarded.

Figure 1 .

 1 Figure 1.15: A sphere rendered by rasterization.

  Figure 1.16: A deferred shading pipeline. It implies two passes of rendering. The rst one to draw geometries in each pixel. The second one to shade each pixel.

Figure 2

 2 Figure 2.1: Examples of interactive renderings produced with 3DCast. Left image shows some transmittance function mapping as described by Delalandre et al. [2011]. Right image: a complex 3D virtual world.

Figure 2 . 2 :

 22 Figure 2.2: Volumetric rendering in 3DCast using the technique introduced by Gautron et al. [2013].

  Figure 2.3: Massive dynamic lighting in 3DCast.

Figure 2 . 4 :

 24 Figure 2.4: Adaptive terrain streaming using 3DCast from the work of Lerbour et al. [2010].

Figure 2 . 5 :

 25 Figure 2.5: Contact visualization in the 3DCast Maya Plugin helps in detecting missing contacts (Right image: red rectangle) and so prevents from oating objects (Left image: nal render with Mental Ray). Middle image shows the standard Maya Viewport.

  4.4. Nevertheless, bidirectional algorithms such as vertex connection and merging by Georgiev et al. [2012a], presented 32 Arthur Dufay 2. Proposed Path Tracing Architecture in 3DCast

Figure 2

 2 Figure 2.6: A diuse blue Phong BRDF left, gold glossy Phong BRDF right.

Figure 2

 2 Figure 2.7: A mix of blue diuse and specular BRDF. With k d = 0.5 and k s = 0.5.

Figure 2 .

 2 Figure 2.10: A statue (Hebemissin model) in a box lit by an area light on the ceiling, casting soft shadows.

  vec3 linearToneMapping ( in vec3 color , in float gammaFactor , in float maxDisplayLuminance ) { // compute luminance of pixel float lum = luminance ( color ) ; // apply gamma correction float gammaCorrectedLuminance = pow ( lum , gammaFactor ) ; vec3 result = vec3 ( gammaCorrectedLuminance / lum ) * color ; // rescale luminance return result / maxDisplayLuminance ; } To compute the luminance of a HDR pixel we use the following equation from Reinhard et al. [2008], that works for a color c in RGB space: luminance = 0.2126 * c.red + 0.7152 * c.green + 0.0722 * c.blue (2.8) 2.2.4 3DCast Path Tracer -Architecture Overview Addressing the Whole GPU Memory: NV_shader_buer_load extension

Figure 2 .

 2 Figure 2.13: Our splatting technique to preview indirect lighting. Following the blue arrow order, rst image: no indirect lighting, second image: only one pixel computed in each tile of 16×16 and splat to all the other pixels in the tiles, third image one pixel in each tile of 8×8 computed, and last image all the pixels computed at least once.

Figure 2 .

 2 Figure 2.14: The six test scenes for the rst frame computation benchmark. From left to right, top to bottom: Siebnik (80K triangles), Interior (85K triangles), Dragon (870K triangles), Museum (1.5M triangles), Cathedral (1.01M triangles) and Hairball (2.88M triangles).

Figure 3

 3 Figure 3.1: Rising of GPU computing power vs CPU power. Courtesy of NVIDIA.

  in one clock cycle. They also have a dedicated mul-add and 48 Arthur Dufay 3. Kernel Implementation of Path Tracing on GPU

Figure 3

 3 Figure 3.2: A view of the NVIDIA GK110 Processor based on Kepler architecture, courtesy of NVIDIA.

•

  Registers 65, i.e., 64KB per SM • L1 cache + Shared Memory 64KB and Read-Only Data Cache 48KB per SM, low latency (10-20 cycles), very high bandwidth 1.5-2.5 TB/s 52 Arthur Dufay 3. Kernel Implementation of Path Tracing on GPU

Figure

  Figure 3.4: NVIDIA general GPU memory architecture

Figure 3 . 5 :

 35 Figure 3.5: Our reverse shadow ray technique launches rays from light sources to surfaces (red arrows).

Figure 3

 3 Figure 3.6: A light sampling buer is constructed, aligned with the image grid. Then, it is split in smaller buers to regroup shadow rays per light source.

LimitationsA

  uniform grid does not adapt to the topology of the scene, thus, leading to an unoptimized memory usage, especially for sparsely distributed geometry See red rectangles in Figure4.1.

Figure 4 . 4 .

 44 Figure 4.1: A uniform grid. The red rectangles show memory waste and unoptimized subdivision leaving to unnecessary traversal.

Figure 4 .

 4 Figure 4.2: Left: 2D example of an octree, called a quadtree. Right: 3D example of an octree.

Figure 4 .

 4 Figure 4.3: Kd-tree examples in 2D (left) and 3D (right).

Figure 4 . 4 :

 44 Figure 4.4: BSP-Tree.

4 .

 4 Analysis for an Adequate Spatial Acceleration Data Structure

Figure 4 .

 4 Figure 4.5: A right side comb shaped BVH.

Figure 4

 4 A path tracing benchmark with and without our Meta-BVH on an NVIDIA GTX Titan X. All scenes were rendered with a max path length of 3 indirect bounces, performances measured on indirect lighting bounces only. The Meta-BVH is protable in pathological cases only, like the Plane Dragons scene.

Figure 4 . 6 :

 46 Figure 4.6: The three test scenes for the Meta-BVH benchmark. From left to right: San Miguel (10M triangles 253 objects), Interior Scene (85K triangles 109 objects), Plane Dragons (872M triangles 1002 objects).

Figure 4 .

 4 Figure 4.8: Left: A goToNextNode() call on a roped BVH, it requires only one node fetch. Right: The new pointer hierarchy (orange curves) of a roped BVH compared to standard

/

  / compute a seed for the RNG based on pixel coordinates void seed ( int x , int y , int screenWidth ) { float seed = 0.0174532 * y * screenWidth + x ; seed *= sin ( seed ) ; return seed ; } vec2 rng ( float & seed ) { // fract returns the fractional part of x . It is calculated as x -floor ( x ) . return fract ( sin ( vec2 ( seed +=0.1 , seed +=0.1) ) * vec2 (43758.5453123 ,22578.1459123) ) ; } However, this does not give really well-distributed samples and is limited in the number of dimensions it can produce (see Figure 5.1).

Figure 5

 5 Figure 5.1: 100 points generated using a pixel-based RNG as described in Algorithm 5.1.

  float wangHash ( uint & seed ) { seed = ( seed ^61) ^( seed >> 16) ; seed *= 9; seed = seed ^( seed >> 4) ; seed *= 0 x27d4eb2d ; seed = seed ^( seed >> 15) ; return ( float ( seed ) ) / 0 xffffffffU ; } An example of a sample set generated with Wang Hash is shown in Figure 5.2. It was rst introduced as a hash function but gives a quite good result as a simple RNG on GPU as presented by[Reed, 2013].

Figure 5

 5 Figure 5.2: Wang Hash, 100 points.

Figure 5

 5 Figure 5.3: Halton sequence, left: dimensions 2 and 3, right: dimensions 7 and 9, 100 points. Observe how sampling quality decreases when getting higher in the dimensions of the sequence, leaving a lot of empty spaces and producing alignments.

Figure 5 . 4 :

 54 Figure 5.4: Sobol sequence, left: dimensions 2 and 3, right: dimensions 7 and 9, 100 points. Observe how the sampling quality decreases when getting higher in the dimensions of the sequence, leaving a lot of empty spaces and producing alignments.

Figure 5

 5 Figure 5.5: Left: coherent path tracing 128spp. Right: path tracing with Cranley Patterson rotation 128spp. Artifacts of coherent path tracing are induced by the re-utilization of the same sampling sequence for each pixel. The structured noise is more prominent on at surfaces.

Figure 5 . 6 :

 56 Figure 5.6: The rotation method to decorrelate samples applies a rotation in disk space (left) to generate new samples in the hemisphere (right).

Figure 5

 5 Figure 5.7: A superposition of a 50 points sample-set of the Halton sequence projected on a disk (black points), and 50 rotated sample-set of this set (blue points). See how the rotation method induces circular patterns.

Figure 5

 5 Figure 5.8: Left: Cranley Patterson applied on two points x 0 and x i of the same input sampling to optain new samples x 0 and x i . As x i falls outside of the domain it gets reprojected using a simple modulo. Right: Cranley Patterson applied on the whole point set.

Figure 5 .

 5 Figure 5.10: Coherent path tracing, which does not use any decorrelation method, presents the best GPU cache coherency and thus the best ray tracing performance.

Figure 5 .

 5 Figure 5.11: Our micro jitter method generates new sample sets (orange and blue) from an input sequence (black) by jittering samples.

Figure 5 .

 5 Figure 5.14: Plot of the star discrepancy of the Halton sequence depending on the number of samples. See how f (x) = 2.5 √ N can t the star discrepancy.

Figure 5 .

 5 Figure 5.15: Plot of the star discrepancy of the Hammersley sequence depending on the number of samples. See how f (x) = 1.5 √ N can t the star discrepancy.

Figure 5 .

 5 Figure 5.16: The three test scenes for our micro jitter performance analysis. From left to right: Interior scene (85K triangles), Dragon scene (870K triangles) and Museum (1.5M triangles).

Figure 5 .

 5 Figure 5.19: Comparison of our µ-jitter method vs the Cranley Patterson rotation on AO computation at dierent spp. Our new decorrelation technique does not exihibit any loss in image quality.

  Figure 1.18) is an interesting research eld. Furthermore, VPL methods are faster to converge than PT methods in at diuse 104 Arthur Dufay Conclusion areas. See Figure 5.21 for a comparison of the two of them.

Figure 5 .

 5 Figure 5.21: Path tracing (50spp) vs 25600 VPLs. Each method uses one bounce of indirect illumination. Path tracing still presents some noise and VPLs display artifacts.

Figure 5 .

 5 Figure 5.22: Combination of VPLs and path tracing for indirect illumination computation. One bounce only. 50 spp for path tracing and 25600 VPLs.

  Figure B.7).

  [2009] we propose to explicitly store geometric information regarding the closest visible triangle in shadow map pixels. To that end we use the compact triangle storage described in Figure B.1.

Figure B. 1 :

 1 Figure B.1: Explicit triangle storage using our compact representation.

  using this representation is that the shadow ray intersection is less expensive compared to shadow ray intersection with full triangle 3D coordinates. It is accomplished using inexpensive early rejection test and simple 2D point in triangle. We start by reconstructing the depth at view sample projection p over the triangle plane using the rst order approximation formula:d Z = p Z + (p x -c x ) d p Z dx + (p y -c y ) d p Z dy(B.1) See Figure B.2 for notations.

Figure B. 2 :

 2 Figure B.2: Triangle notations.

Figure B. 3 :

 3 Figure B.3: Incorrect depth spawned using slope scale minimum depth.

  Figure B.4 case a).

Figure B. 4 :

 4 Figure B.4: Minimum depth nding cases.

  Figure B.4 case b).

  Figure B.4 case c).depth min = p Z + 0.5 d p Z dx + 0.5 d p Z dy (B.3)This corresponds to the same formula used by the slope scale based depth correction.

  Figure B.4 case d).

Figure B. 5 :

 5 Figure B.5: Conservative shadow map generation process.

Figure B. 6 :

 6 Figure B.6: Lookup pixel neighborhood to enforce pixel classication and avoid ray tracing operations.

  the minimum depth stored in the pixel is greater or nearly equal to the depth of the view sample in light space, then the view sample (or pixel) is classied as lit.

Figure B. 7 :

 7 Figure B.7: Our improved classication drastically reduces the number of ray traced shadows (red pixels).

B. 5 .

 5 Advantages of the Invention the depth buer. Thanks to the conservative rasterization we are sure that no other geometry is present between beyond the ray.

Figure

  Figure B.8: GPU ray traced shadow using the minimum depth shadow map.

Table 2

 2

	Company Name	Renderer Name Rendering Algorithms
	Pixar [2017]	Renderman	Path Tracing -Vertex Connec-
			tion and Merging (VCM)
	Solid Angle [2017]	Arnold	Path Tracing
	Chaos Group [2017]	V-Ray	Path Tracing
	NVIDIA [2017b]	Mental Ray	Path Tracing
	MaxwellRender [2017]	Maxwell	Hybrid of Bidirectional Path
			Tracing and Metropolis Light
			Transport

.1: Some of the most used production renderers for VFX with their corresponding rendering algorithms.

Table 2 .

 2 2: Benchmark of several implementations of our path tracer, images of 1280x720

	2. Proposed Path Tracing Architecture in 3DCast
			One-Pass	16-Pass blocks	256-Pass Quad-tree
	Scene	SPP	FPS Time (s) FPS	Time (s) Time Ratio	FPS	Time (s) Time Ratio
	Siebnik	50	0.20	254.2	2.91	275.1 1.08	35.96	356.0 1.40
	Interior	100 0.31	319.5	4.56	351.1 1.10	46.49	550.7 1.72
	Dragon	100 1.47	67.93	18.62	85.91 1.26	128.26	199.6 2.94
	Cathedral 100 0.76	131.4	10.13	157.9 1.20	65.37	391.6 2.98
	Museum	50	0.14	347.9	2.11	378.9 1.09	26.82	477.3 1.37
	Hairball	50	0.22	227.3	2.93	272.6 1.20	18.31	699 3.08

Table 3 .

 3 1: Terminology of the CUDA programming language.

	SIMT	Single Instruction Multiple Threads, see Flynn
		[1972] taxonomy
	Workgroup/Block Logical structure that groups threads together
	Warp	Group of 32 threads executed in an SIMT manner
	SM	Streaming Multiprocessor, physical structure that
		executes one or more Workgroups
	SMX	New generation SM
	GPC	Graphics Processing Cluster, physical structure
		that groups SMs together
	Grid	Whole set of threads that execute the kernel, con-
		tains one or more Workgroups
	Kernel	Set of instructions: code to be executed on the GPU
		by the Workgroups
	SFU	Special function unit, executes specic functions
		(cosinus, log, ...)
	ALU	Integer Arithmetic Logic Unit
	FPU	Floating Point Unit

  4. Analysis for an Adequate Spatial Acceleration Data Structure

Table 5 .

 5 1: Results of our RMSE measurements on ambient occlusion images generated with Cranley Patterson rotation method (CP) and our micro jitter, with Halton sequence and Hammersley sequence. Red numbers highlight dierences between CP and ours.

	RMSE -AO -Halton Sequence	RMSE -AO -Hammersley Sequence
	spp	CP	Micro-jitter	spp	CP	Micro-jitter
	8	0.0994	0.0990	8	0.0950	0.0950
	16 0.0621	0.0624	16 0.0567	0.0572
	32 0.0385	0.0381	32 0.0337	0.0336
	64 0.0236	0.0239	64 0.0200	0.0196
	128 0.0148	0.0148	128 0.0119	0.0116
	256 0.0086	0.0088	256 0.0070	0.0070
	512 0.0058	0.0058	512 0.0042	0.0042

For instance the latest professional NVIDIA card, the P6000 has

24GB of GRAM, whereas CPU can handle up to 512GB of RAM (8x64GB).

Theoretical Background

Proposed Path Tracing Architecture in 3DCast

Arthur Dufay 2. Proposed Path Tracing Architecture in 3DCast

5.4. A GPU Cache Friendly Decorrelation Technique -Micro Jitter

Remerciements

To better understand our implementation choices we need to understand how a GPU works. In this chapter, we present more in detail the GPU architecture and how to implement a path tracer on the GPU.

In Section 3.1, we start by presenting the architecture of a modern GPU and give the basic concepts of general-purpose processing on graphics processing units (GPGPU). Even though our implementation relies on OpenGL Compute Shaders, we base our presentation on the CUDA specication since they share the same philosophy and restrictions on thread assignment and memory limitations.

Then, to review the impact on performance due to implementation choices, we implemented two dierent path tracers that are presented in Sections 3. 4.2 and 3.4.3. Especially, we present in these sections how to split the GPU kernel, and expose a benchmark of these two implementations on several GPUs.

Finally we describe the Fast visibility test using reversed shadow rays patent we submitted, that helps to accelerate shadow-ray queries in Section 3.6.

and BVH

The performance dierences between BVH and kd-trees have been well studied by [START_REF] Vinkler | Performance comparison of bounding volume hierarchies and kd-trees for gpu ray tracing[END_REF]. They studied the ray tracing performance of these two SADS on the GPU over 16 dierent scenes with a wide diversity of geometry complexity. This study proves that BVH performs better on scenes with small to medium sizes (80K triangles to 7M triangles). For larger scenes, kd-trees perform better. This is mostly due to the duplication of triangles that appears inherently in the kd-trees and that does not seem to be an advantage on small to medium scenes. On larger scenes, even though kd-trees duplicate a lot of references to triangles, they allow to discard a lot of 3D space and thus give better performances than BVHs.

However, in our context, we thought that the BVH is the most appropriated data structure, for its small memory consumption and capability to handle relatively complex scenes. At this time, GPU memory is still smaller than CPU memory and we need to be as compact as possible. For instance the latest professional NVIDIA card, the P6000 has 24GB of GRAM, whereas CPU can handle up to 512GB of RAM (8 x 64GB). Since we chose to implement the BVH in our path tracer, we present in more details how we did to intersect it on the GPU. 5.1). Our technique did not show any loss in image quality (cf. Figure 5.19).

Application to Screen Space Sampling

Even though our µ-jitter was designed for decorrelating path tracing samples, we found an application for faster sampling in screen space techniques. Our method can be applied to screen space ambient occlusion (SSAO) computation for instance. Indeed, SSAO computation requires fetching several samples in a 2D texture around a pixel to average depth values stored in a Z-Buer. These samples lie in a xed size kernel and have a randomly chosen position. In most cases the random position of the pixels samples are precomputed, a Poisson sampling is often used. To apply our method in this case, one can reuse the same sample set across all pixels, and alter it with a constrained jitter vector. The sample set corresponds to the 2D position of the fetched pixels.

We tested it and got faster computation of SSAO. On average, we saw an acceleration factor of 1.47.

Software Tools

During this thesis, two tools outside of 3DCast were developed They are presented here.

A.1 HDR Viewer

We can dump frames in raw HDR (RGBA 32 bits per channel) from 3DCast. This snapshots can then be opened in our HDR Viewer tool (cf.

B.2.1 Hybrid GPU Pipeline for Alias Free Shadows

The closest invention related to ours is the hybrid GPU rendering pipeline for Alias-Free Shadows proposed by [START_REF] Hertel | A hybrid gpu rendering pipeline for alias-free hard shadows[END_REF]. They combine rasterization of a shadow map with a GPU ray tracer to accelerate the rendering of ray traced shadows. To that end they consider the conservative rasterization of a shadow map in which they store depth, triangle ID and a ag indicating if the triangle fully covers or not the shadow map pixel. Using a simple depth test and checking the triangle coverage ag, they are able to quickly classify pixels