Index of /rozprawy2/10288

Pełen tekst

(1)AGH University of Science and Technology Faculty of Electrical Engineering, Automatics, Computer Science and Electronics Institute of Computer Science. Dissertation for the degree of Doctor of Philosophy in Computer Science. Design and Implementation of Parallel Light Transport Algorithms based on quasi-Monte Carlo Ray Tracing. mgr in˙z. Michal Radziszewski. supervisor: dr hab. in˙z. Krzysztof Boryczko, prof. n. AGH. Krak´ow, June 2010.

(2) Abstract Photorealistic rendering is a part of computer graphics, which concentrates on creating images and animations based on 3D models. Its goal is creation of pictures that are indistinguishable from real world scenes. This work is dedicated to a particular class of photorealistic rendering algorithms – quasi-Monte Carlo ray tracing based global illumination. Global illumination is a very useful concept for creating realistically lit images of artificial 3D scenes. Using automatic and correct computation of vast diversity of optical phenomena, it enables creating a rendering software, which allows specification of what is to be rendered instead of detailed description of how to render a given scene. Its current applications range from many CAD systems to special effects in movies. In future, when computers become sufficiently powerful, real time global illumination may be the best choice for computer games and virtual reality. Currently, only Monte Carlo and quasi-Monte Carlo ray tracing based algorithms are general enough to support full global illumination. Unfortunately, they are very slow compared to other techniques, e.g. hardware accelerated rasterization. The main purpose of this thesis is an improvement of efficiency of physically correct rendering. The thesis concentrates on enhancement of robustness of rendering algorithms, as well as parallel realization of them. These two elements together can substantially increase global illumination applicability, and are a step towards the ultimate goal of being able to run true global illumination in real time.. i.

(3) Acknowledgements I am deeply indebted to my supervisor, Prof. Krzysztof Boryczko. Without his inspiration, advice and encouragement this work would not have been completed. I would like to sincerely thank Dr Witold Alda from the AGH University of Science and Technology. During work on this dissertation we have spent many hours in conversations. The expertise I gathered working with him contributed substantially to the dissertation. I would like to acknowledge the financial support from the AGH University of Science and Technology – a two one year scholarships for PhD students. Finally, I also would like to acknowledge ZPORR (Zintegrowany Program Operacyjny Rozwoju Regionalnego) for scholarship “Malopolskie Stypendium Doktoranckie”, co-founded by European Union.. ii.

(4) Contents Abstract. i. Acknowledgements. ii. 1 Introduction 1.1 Global Illumination and its Applications . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Thesis Purpose and Original Contributions . . . . . . . . . . . . . . . . . . . . . . 1.3 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1 1 2 3. 2 Light Transport Theory 2.1 Geometric Optics . . . . . . . . . . . . . . 2.1.1 Assumptions . . . . . . . . . . . . 2.1.2 Radiometric Quantities . . . . . . 2.2 Light Transport Equation . . . . . . . . . 2.2.1 Surface Only Scattering . . . . . . 2.2.2 Volumetric Scattering Extension . 2.2.3 Properties of Scattering Functions 2.2.4 Analytic Solutions . . . . . . . . . 2.2.5 Simplifications . . . . . . . . . . . 2.3 Image Formation . . . . . . . . . . . . . . 2.3.1 Importance . . . . . . . . . . . . . 2.3.2 Integral Formulation . . . . . . . . 2.3.3 Image Function . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. 4 4 4 5 7 7 7 9 10 10 11 11 12 13. 3 Monte Carlo Methods 3.1 Monte Carlo Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Statistical Concepts . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Estimators of Integrals . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Biased and Unbiased Methods . . . . . . . . . . . . . . . . . . . 3.2 Variance Reduction Techniques . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Importance and Multiple Importance Sampling . . . . . . . . . . 3.2.2 Russian Roulette and Splitting . . . . . . . . . . . . . . . . . . . 3.2.3 Uniform Sample Placement . . . . . . . . . . . . . . . . . . . . . 3.3 Quasi-Monte Carlo Integration . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Desired Properties and Quality of Sample Sequences . . . . . . . 3.3.2 Low Discrepancy Sequences . . . . . . . . . . . . . . . . . . . . . 3.3.3 Randomized Quasi-Monte Carlo Sampling . . . . . . . . . . . . . 3.3.4 Comparison of Monte Carlo and Quasi-Monte Carlo Integration . 3.3.5 Quasi-Monte Carlo Limitations . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. . . . . . . . . . . . . . .. 14 14 14 16 16 17 17 18 18 19 20 21 23 23 24. 4 Light Transport Algorithms 4.1 Ray Tracing vs. Other Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 View Dependent vs. View Independent Algorithms . . . . . . . . . . . . . . 4.1.2 Ray Tracing Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 26 27 27 28. iii. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . ..

(5) CONTENTS. 4.2. 4.3. 4.4. 4.5. iv. 4.1.3 Hardware Accelerated Rasterization . . . . . . . . . . . . 4.1.4 Radiosity Algorithms . . . . . . . . . . . . . . . . . . . . Light Transport Paths . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Classification of Paths . . . . . . . . . . . . . . . . . . . . 4.2.2 Construction of Paths . . . . . . . . . . . . . . . . . . . . 4.2.3 Local Path Sampling Limitation . . . . . . . . . . . . . . Full Spectral Rendering . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Necessity of Full Spectrum . . . . . . . . . . . . . . . . . 4.3.2 Representing Full Spectra . . . . . . . . . . . . . . . . . . 4.3.3 Efficient Sampling of Spectra . . . . . . . . . . . . . . . . 4.3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . Analysis of Selected Light Transport Algorithms . . . . . . . . . 4.4.1 Path Tracing . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Bidirectional Path Tracing . . . . . . . . . . . . . . . . . 4.4.3 Metropolis Light Transport . . . . . . . . . . . . . . . . . 4.4.4 Irradiance and Radiance Caching . . . . . . . . . . . . . . 4.4.5 Photon Mapping . . . . . . . . . . . . . . . . . . . . . . . Combined Light Transport Algorithm . . . . . . . . . . . . . . . 4.5.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.2 Merging of an Unbiased Algorithm with Photon Mapping 4.5.3 Results and Conclusion . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. 28 29 29 30 30 31 32 33 35 36 39 42 42 46 49 51 52 57 57 58 59. 5 Parallel Rendering 5.1 Stream Processing . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Stream Processing Basics . . . . . . . . . . . . . . . 5.1.2 Extended Stream Machines with Cache . . . . . . . 5.1.3 Stream Monte Carlo Integration . . . . . . . . . . . 5.2 Parallel Ray Tracing . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Algorithm Initialization and Scene Description . . . 5.2.2 Frame Buffer as an Output Stream . . . . . . . . . . 5.2.3 Multipass Rendering . . . . . . . . . . . . . . . . . . 5.2.4 Ray Tracing on an Extended Stream Machine . . . . 5.2.5 Results and Conclusions . . . . . . . . . . . . . . . . 5.3 Choice of Optimal Hardware . . . . . . . . . . . . . . . . . 5.3.1 Shared Memory vs. Clusters of Individual Machines 5.3.2 Multiprocessor Machines vs. Graphics Processors . . 5.3.3 Future-proof Choice . . . . . . . . . . . . . . . . . . 5.4 Interactive Visualization of Ray Tracing Results . . . . . . . 5.4.1 Required Server Output . . . . . . . . . . . . . . . . 5.4.2 Client and Server Algorithms . . . . . . . . . . . . . 5.4.3 MIP-mapping Issues . . . . . . . . . . . . . . . . . . 5.4.4 Results and Discussion . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . .. 61 62 62 62 64 65 65 65 66 66 68 69 70 70 71 71 72 74 78 79. 6 Rendering Software Design and Implementation 6.1 Core Functionality Interface . . . . . . . . . . . . . 6.1.1 Quasi-Monte Carlo Sampling . . . . . . . . 6.1.2 Ray Intersection Computation . . . . . . . 6.1.3 Spectra and Colors . . . . . . . . . . . . . . 6.1.4 Extension Support . . . . . . . . . . . . . . 6.2 Procedural Texturing Language . . . . . . . . . . . 6.2.1 Functional Languages . . . . . . . . . . . . 6.2.2 Syntax and Semantic . . . . . . . . . . . . . 6.2.3 Execution Model and Virtual Machine API 6.2.4 Results and Conclusion . . . . . . . . . . . 6.3 New Glossy Reflection Models . . . . . . . . . . . . 6.3.1 Related Work . . . . . . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. 82 82 83 85 86 89 90 91 92 93 95 96 97. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . ..

(6) CONTENTS. v. 6.3.2 6.3.3 6.3.4. Properties of Reflection Functions . . . . . . . . . . . . . . . . . . . . . . . 98 Derivation of Reflection Function . . . . . . . . . . . . . . . . . . . . . . . . 99 Results and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101. 7 Results 104 7.1 Image Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.2 Full Spectral Rendering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.3 Comparison of Rendering Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 105 8 Conclusion 109 8.1 Contributions Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 8.2 Final Thoughts and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Bibliography. 112. Index. 121.

(7) List of Symbols Φ E L W λ fr ft fs fp N ω ωi ωo θ x, y x\ 2 − x1 x, x[k] Λ Ω Ω+ A V X µ σ(ω) σ ⊥ (ω) A(x) V (x) σa σs σe pdf cdf ξ δ. radiant flux irradiance radiance importance wavelength bidirectional reflection distribution function (BRDF) bidirectional transmission distribution function (BTDF) bidirectional scattering distribution function (BSDF) phase function normal direction unit directional vector incident ray direction outgoing ray direction angle between ω and N points – either on a surface A or in a volume V normalized vector pointing from x1 to x2 , x1 6= x2 a light transport path, path with k segments, k + 1 vertexes space of all visible wavelengths, Λ = [λmin , λmax ] space of all unit directional vectors ω space of all unit directional vectors such that ω ◦ N ≥ 0 space of all points on scene surfaces space of all points in scene volume, without surfaces space of all light transport paths arbitrary measure solid angle measure projected solid angle measure area measure volumetric measure absorption coefficient scattering coefficient extinction coefficient probability density function cumulative distribution function canonical uniform random variable Dirac delta distribution. vi.

(8) List of Figures 2.1 2.2. Radiance definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results of simplifications of light transport equation . . . . . . . . . . . . . . . . .. 6 11. 3.1 3.2 3.3 3.4 3.5. Uniform sample placement . . . . . . . . . . . . . . Quasi-Monte Carlo sampling patterns . . . . . . . Comparison of Monte Carlo and quasi-Monte Carlo Undesired correlations between QMC sequences . . Rendering with erroneous QMC sampling . . . . .. 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 4.18 4.19 4.20 4.21. An example light path. . . . . . . . . . . . . . . . . . . . . Extended light path notation . . . . . . . . . . . . . . . . Difficult light path . . . . . . . . . . . . . . . . . . . . . . Local path sampling limitation . . . . . . . . . . . . . . . Full spectral and RGB reflection . . . . . . . . . . . . . . Full spectral and RGB refrefraction on prism . . . . . . . Different methods for wavelength dependent scattering . . Selection of optimal number of spectral samples . . . . . . Various methods of sampling spectra . . . . . . . . . . . . Imperfect refraction with dispersion . . . . . . . . . . . . Analysis of behaviour of spectral sampling . . . . . . . . . A path generated by Path Tracing algorithm . . . . . . . Results of Path Tracing . . . . . . . . . . . . . . . . . . . Simplified Path Tracing . . . . . . . . . . . . . . . . . . . Batch of paths generated with BDPT. . . . . . . . . . . . Batch of paths generated with optimized BDPT. . . . . . Comparison of BDPT and Photon Mapping . . . . . . . . One pass versus two pass Photon Mapping . . . . . . . . . Quickly generated images with one pass Photon Mapping Light transport algorithms comparison . . . . . . . . . . . Light transport algorithms comparison . . . . . . . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. 20 23 24 25 25. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. 29 31 31 32 34 34 38 38 39 40 41 43 45 45 46 48 55 56 57 59 60. 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10. Stream machine basic architecture. . . . . . . . . . . . . . . . . . . . Extended stream machine architecture. . . . . . . . . . . . . . . . . . Test scenes for parallel rendering . . . . . . . . . . . . . . . . . . . . Parallel rendering run times . . . . . . . . . . . . . . . . . . . . . . . Main loop of visualization client process. . . . . . . . . . . . . . . . . Glare effect applied as a postprocess on the visualization client. . . . Comparison of MIP-mapping and custom filtering based blur quality. Results of Interactive Path Tracing. . . . . . . . . . . . . . . . . . . Results of Interactive Photon Mapping. . . . . . . . . . . . . . . . . Noise reduction based on variance analysis of Path Tracing. . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . .. 62 63 68 68 75 77 79 79 80 80. 6.1 6.2. Ray-primitives interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Semi-transparent surfaces intersection optimization . . . . . . . . . . . . . . . . . .. 86 86. vii. . . . . . . . . . . . . . . . . . . . . integration error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . ..

(9) LIST OF FIGURES. viii. 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14. Comparison between different gamut mapping techniques Textured and untextured models . . . . . . . . . . . . . . Sample procedural texture. . . . . . . . . . . . . . . . . . Images generated using noise primitive . . . . . . . . . . . Procedurally defined materials . . . . . . . . . . . . . . . Mandelbrot and Julia Fractals . . . . . . . . . . . . . . . . Comparison of different glossy BRDFs with little gloss . . Latitudal scattering only . . . . . . . . . . . . . . . . . . . Longitudal scattering only . . . . . . . . . . . . . . . . . . Product of latitudal and longitudal scattering . . . . . . . Scattering with perpendicular and grazing illumination . . Complex dragon model rendered with our glossy material. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. 89 90 91 95 96 96 102 102 102 102 103 103. 7.1 7.2 7.3. Comparison of spectral rendering algorithms . . . . . . . . . . . . . . . . . . . . . . 105 Full spectral rendering of a scene with imperfect refraction . . . . . . . . . . . . . . 107 Rendering of indirectly visible caustics . . . . . . . . . . . . . . . . . . . . . . . . . 108.

(10) List of Tables 4.1. Numerical error of spectral sampling . . . . . . . . . . . . . . . . . . . . . . . . . .. 40. 6.1 6.2 6.3. Spectral functions for RGB colors . . . . . . . . . . . . . . . . . . . . . . . . . . . . Texturing language grammar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Notation used in BRDF derivation. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 88 94 99. 7.1 7.2. Convergence of spectral sampling techniques . . . . . . . . . . . . . . . . . . . . . . 105 Convergence of selected rendering algorithms . . . . . . . . . . . . . . . . . . . . . 106. ix.

(11) List of Algorithms 4.1 4.2 5.1 5.2 5.3 5.4 6.1 6.2. Optimized Metropolis sampling . . . . . . . . . . . . . . . . . Construction of photon maps . . . . . . . . . . . . . . . . . . Kernel for Monte Carlo integration . . . . . . . . . . . . . . . Single- and multipass rendering on extended stream machine. Rasterization of point samples by Visualization client . . . . Visualization client repaint processing . . . . . . . . . . . . . Gamut mapping by desaturation . . . . . . . . . . . . . . . . Sample usage of procedural texturing . . . . . . . . . . . . .. x. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. 51 53 64 66 76 77 89 95.

(12) Chapter 1. Introduction Photorealistic rendering is a part of computer graphics, which concentrates on creating still images and animations based on 3D models. Its goal is creation of pictures that are indistinguishable from real world scenes. This work is dedicated to the particular class of photorealistic rendering algorithms – ray tracing based global illumination. In the rest of this chapter, a brief description of global illumination is presented, followed by the outline of the most interesting original contributions, and finally, the thesis organization.. 1.1. Global Illumination and its Applications. The global illumination term is used to name the two distinct types of capabilities of rendering algorithms – there are two commonly used and substantially different definitions of it. According to the first definition, global illumination effects are just opposed to local illumination effects, where only direct lighting is accounted for while rendering any part of the scene. Thus, any phenomenon which depends on knowledge of other scene parts, while rendering a given primitive, is a global illumination effect. Such effects are, for example, shadow casting and environment mapping. On the other hand, according to the second definition, any rendering algorithm capable of global illumination must be able to simulate all possible interreflections of light between scene primitives. Since the second definition is much more useful and precise, it is used through the rest of the thesis. Global illumination is a very useful concept for creating realistically lit images of artificial 3D scenes. Computing automatically and correctly vast diversity of optical phenomena, it creates a solid basis for rendering software, which allows specification of what is to be rendered instead of detailed description of how to render a given scene. Its current applications range from many CAD systems to special effects in movies. Global illumination algorithms are responsible for the rendering process. Within the framework presented later they are easy to implement, however design of an effective, robust and physically correct global illumination algorithm is a difficult and still not fully solved task. By physical correctness we understand complete support of geometric optics based phenomena. Currently, only ray tracing based algorithms are general enough to render all geometric optic effects. Unfortunately, the price to pay for such an automatization is that the evaluation of global illumination is slow compared to other techniques, e.g. hardware accelerated rasterization. However, it is widely believed that, when computers become fast enough, global illumination is likely to replace other, less physically accurate techniques. This statement is often supported by the fact that a similar breakthrough is already happening in modeling domain, where physical models compete with more traditional approaches with good effect. For example, nowadays it is possible to model a cloth animation either as a mesh with a given mass, stiffness, etc. and let the model to calculate its appearance over time, or as a sequence of keyframes, each frame laboriously. 1.

(13) CHAPTER 1. INTRODUCTION. 2. created by animator. Currently both models can be run in real time, while the simulation was not plausible a few years ago.. 1.2. Thesis Purpose and Original Contributions. The main purpose of this thesis is to improve an efficiency of physically correct rendering. The idea is to provide working software with rigorous theoretical basis. Our assumption is to avoid two extremes – overly theoretical work, without care about algorithms implementation and algorithms created by method of trials and errors, designed to produce good enough visual effects, which sometimes work and sometimes do not. The latter approach is, unfortunately, surprisingly common in real time graphics, despite it does not provide any information about algorithm correctness. Moreover, accounting for hardware development trends is very important for us, since this can significantly affect performance of algorithms. The thesis concentrates on enhancement of robustness of rendering algorithms, as well as their parallel realization. These two elements together can substantially increase global illumination applicability, and are a step towards the ultimate goal of being able to run true global illumination in real time. In effort to realize it, this dissertation provides several original contributions in the field of computer graphics. This section enumerates most important and interesting of them, listed in order in which they appear in the thesis. Light transport theory. Typically light transport theory assumes creation of image from pixels by convolution of luminance with filter function associated with each pixel. There is nothing incorrect with this concept, but it limits generality of rendering algorithms. Instead, we represent image as a 3D function defined over [0, 1]2 × λ, where unit square represents image film surface and λ is wavelength of light. This generalization allows using much more sophisticated post-processing techniques, but it makes invalid any rendering algorithm dependent on pixel abstraction, however. This approach is explained in details in Section 2.3. Full spectral rendering. Despite being proven incorrect, a lot of global illumination algorithms are designed to use an RGB color model. Only few most advanced approaches attempt to accurately simulate visually pleasing full spectral phenomena. We have designed and implemented an improved full spectral support, based on Multiple Importance Sampling. This new technique is much more efficient at simulating non-idealized wavelength dependent phenomena, and fits elegantly into Monte Carlo and quasi-Monte Carlo ray tracing algorithms. The novel spectral sampling technique is defined in Section 4.3. Light transport algorithms. First, we have adapted Photon Mapping to be one pass technique. The modified algorithm starts storing just few photons, and later, the photon map size is increased. This aspect allows rendering images with progressive quality improvement, a feature impossible to obtain in the original, two pass variants of this technique. Surprisingly, dynamic photon map comes at a little performance penalty for typical rendered scenes. Second, we have provided a detailed potential error prediction technique for some most significant ray tracing algorithms. This feature is then used in the presented new rendering algorithm, which tries to select the most appropriate method to render any part of an image. Both algorithms are presented in Sections 4.4 and 4.5. Parallel rendering. We provide an extension to the stream processor model to support readwrite memory. The extension guarantees coherency of all pieces of written data, but the order of different reading and writing operations is not preserved. Therefore, the correctness of any algorithm must not depend on content of this kind of memory, but the algorithm may use it to accelerate its operations. This mechanism is used in parallel implementation of one pass version of photon mapping, as well as in our new rendering algorithm. The extended stream machine is described in Chapter 5. Moreover, we have designed and implemented an interactive viewer of.

(14) CHAPTER 1. INTRODUCTION. 3. ray tracing results, based on processing power of GPUs. The viewer works in parallel with CPU based renderer, which generates new data while previous is displayed. This concept is explained in Section 5.4. Sampling oriented interface. We have designed an interface between 3D objects, cameras and rendering algorithms based entirely on sampling. This interface provides clear abstraction layer, general enough to express majority of ray tracing algorithms. It is designed to simplify the implementation of bidirectional light transport methods. Furthermore, the interface provides support for spectral rendering and carefully designed quasi-Monte Carlo sequence generation infrastructure. Additionally, we have developed a technique for storing 2D surfaces and 3D participating media in the same ray intersection acceleration structure, using the sampling interface. If a rendered scene contains roughly similar number of different 2D and 3D entities, the improved algorithm is nearly twice as fast as algorithm using two separate structures. The design of this interface is explained in Section 6.1. Materials. We have provided a shading language optimized for ray-tracing. The new concept is based on usage of a functional language for this purpose. The language is computationally complete and enables easy creation of complex material scripts. The script style resembles much more mathematical notation than classic imperative programming languages. The language is presented in Section 6.2. Moreover, we have developed a new glossy reflection model, which is together symmetric and energy preserving. Derivation of its formulae is given in Section 6.3.. 1.3. Thesis Organization. The second and third chapters provide a theoretical basis used in the rest of the thesis. The second chapter presents a brief introduction to the light transport theory under the assumption of geometric optics applicability. It explains how illumination in a 3D scene can be described mathematically as an integral equation (so called Light Transport Equation), and how to use its solution to create images. The third chapter describes Monte Carlo integration as a general purpose numerical technique, giving details on selected approaches to improve its efficiency. The fourth chapter shows how to apply the Monte Carlo integration to solve Light Transport Equation, which leads to variety of so called non-deterministic ray tracing algorithms. The main point of this chapter is, however, the analysis of strong and weak points of major existing light transport algorithms, as well as proposal of new algorithm, designed to efficiently cope with many of their issues. Finally, this chapter expresses how efficiently incorporate full spectral rendering into presented methods. The fifth chapter illustrate the potential of parallel rendering. Furthermore, it gives the idea of how to express ray tracing as an algorithm dedicated to slightly extended streaming machine, and describes variety of hardware, as potential candidates for implementation of streaming machine. In this chapter there is also presented an interactive previewer of ray tracing results. The sixth chapter presents the design of rendering software for developing and evaluating light transport algorithms, and some most interesting implementation details. The seventh chapter provides the most important results. Then, it contains detailed comparison of efficiency of all presented light transport algorithms. The last chapter summarizes the original contributions, and finally, it presents some ideas of future work dedicated to rendering..

(15) Chapter 2. Light Transport Theory The light transport theory provides a theoretical basis for an image creation process. The theory used in this thesis is based on the assumption of geometric optics applicability. It describes equilibrium radiance distribution over an entire scene, and additionally, a method for the calculated radiance to form a final image. In computer graphics, light transport was first described formally by Kajiya [Kajiya, 1986]. More recent works, which cover this area, are [Veach, 1997] and [Pharr & Humphreys, 2004]. Arvo’s course [Arvo, 1993] and Pauly’s thesis [Pauly, 1999] provide a detailed description of light transport in volumes. This chapter starts with explanation of assumptions of geometric optics. Next, it presents the equation describing radiance distribution in both forms – for surfaces placed in vacuum only, and for volumetric effects as well. Because all these equations are commonly known, they are discussed briefly. Finally, it is shown how computed radiance is used to form a final image. The presented image formation theory is a modified approach, designed to remove assumption that image is built from pixels.. 2.1. Geometric Optics. There are a few physical theories describing the behaviour of light. Most important of them are: geometric optics, wave optics and quantum optics. In simplification, geometric and wave optics explain transport of light and quantum optics describes its interaction with matter. Each of these theories predicts selected real-world phenomena with certain degree of accuracy, and likely none of them is completely correct. Choice of a physical theory for simulation is a tradeoff between the desired accuracy of a solution and a computational cost. An interesting fact is that, for computer graphics needs, the simplest theory, geometric optics, typically is sufficient to provide high degree of realism. Geometric optics based rendering is perfectly fine at capturing phenomena such as soft shadows, indirect lighting, dispersion, etc. However, despite its physical simplicity, very few rendering systems provide full support of geometric optics based rendering, and any application, which attempts to do so, is far too slow to be used in real time.. 2.1.1. Assumptions. The theoretical model of geometric optics is based on the following simplifying assumptions, which significantly improve the computation efficiency and still allow simulation of majority of commonly seen phenomena: • the number of photons is huge while the photon energies are extremely small – any distribution of photons may be treated as a continuous value; • photons do not interact with each other, thus effects such as interference cannot be simulated; 4.

(16) CHAPTER 2. LIGHT TRANSPORT THEORY. 5. • photon collisions with surfaces and particles of non transparent volumes (e.g. fog or dust) are elastic, which means that photons cannot change wavelength during scattering; • diffraction, continuously varying refractive indexes and all other phenomena which could affect movement of photons are neglected, so between collisions photons travel along straight lines; • speed of photons is infinitely large, the scene is assumed to be always in equilibrium state; • optical properties of materials do not depend on illumination power, therefore illumination is linear, i.e. it can be computed independently for each light source and summed to form final result. Using all these assumptions light transport can be described by means of radiometric quantities, such as flux or radiance. Some phenomena, which require solving wave equation of light transport, like diffraction or interference, cannot be simulated using these quantities at all. However, some other selected non-geometric effects can be easily simulated, by simple extension of these quantities. For example, spectral radiance is used to simulate wavelength dependent effects. In a similar way radiance can be extended to support polarization. Moreover, scattering can be extended to support fluorescence. If spectral radiance is represented as a vector, then reflectivity can be represented as a matrix and scattering event as a matrix-vector multiplication. The simplified case of elastic photon scattering is therefore represented with diagonal matrixes. Implemented software extends geometric optics to support spectral radiance only.. 2.1.2. Radiometric Quantities. Under the assumption of applicability of geometric optics it is enough to use a few radiometric quantities to fully describe light transport in any 3D scene. Each quantity is defined by measuring the distribution of radiation energy with respect to some parameters. Any of these quantities can be defined in standard and spectral (denoted by dependence on λ) version. Radiant Flux Radiant flux is defined as radiated energy per time: Φ(t) =. dQ(t) , dt. Φλ (t, λ) =. d2 Q(t, λ) . dtdλ. (2.1). Radiant flux is measured in Watts. This quantity is useful for description of total emission from a 3D object or a single light source. Irradiance Irradiance in a point x is defined as radiant flux per area: E(x) =. dΦ(x) , dA(x). Eλ (x, λ) =. d2 Φλ (x, λ) . dA(x)dλ. (2.2). Irradiance is measured in Watts per square meter. It is used to describe how strong the illumination on a given surface is..

(17) CHAPTER 2. LIGHT TRANSPORT THEORY. 6. Radiance Radiance is considered to be the most basic quantity in radiometry, and is defined as a power per area per solid angle (see Figure 2.1): L(x, ω) =. d2 Φ(x, ω) , dA⊥ (x)dσ(ω). Lλ (x, ω, λ) =. d3 Φλ (x, ω, λ) , dA⊥ (x)dσ(ω)dλ. (2.3). where dA⊥ (x) is a projected area measure. Radiance is measured in watts per square meter per steradian. It may be rewritten to the more convenient expression: L(x, ω) =. d2 Φ(x, ω) d2 Φ(x, ω) = , |ω ◦ N (x)|dAdσ(ω) dA(x)dσ ⊥ (ω). (2.4). which uses standard area measure on surfaces. Light transport equations are based on radiance, which has a useful property – it is constant when light travels along straight lines in vacuum. To render an image it is enough to know the radiance on camera lens, however some techniques try to compute radiance everywhere. During scattering of photons it is important to distinguish between incident and outgoing radiances. These quantities (defined on the same x) are often marked as Li and Lo .. N. w. q. dw. dAcosq dA Figure 2.1: Radiance of conical bundle of rays is defined as derivative of bundle power with respect to angular divergence dω of the cone and perpendicular area dA of its base. The radiance can be measured on surface non-perpendicular to bundle direction, by projecting the measured area.. Radiant Intensity Radiant intensity is defined as power per unit angle: I(ω) =. dΦ(ω) , dσ(ω). Iλ (ω, λ) =. d2 Φλ (ω, λ) . dσ(ω)dλ. (2.5). Radiant intensity is measured in watts per steradian. It is used to describe how strong the illumination in a particular direction is. In strictly physically based systems, the radiant intensity from any particular point x is always zero, however the I(x, ω) quantity is useful in description of emission from, commonly used in modeling, point-based light sources. Volume Emittance Volume emittance is similar to radiance, but is defined with respect to volume, instead of surface: Lv (x, ω) =. d2 Φ(x, ω) , dV (x)dσ(ω). Lv,λ (x, ω, λ) =. d3 Φλ (x, ω, λ) , dV (x)dσ(ω)dλ. (2.6). where V (x) means volume measure. Volume emittance is measured in Watts per cubic meter per steradian. This quantity is used to describe volumetric phenomena, for example emission from correctly modeled fire..

(18) CHAPTER 2. LIGHT TRANSPORT THEORY. 2.2. 7. Light Transport Equation. All the rendering algorithms have to solve the light transport problem. In 1986 J. Kajiya [Kajiya, 1986] first noticed that his theoretically correct (under the assumption of geometric optics applicability) approach solves equation described in his paper, and all currently available algorithms make simplifications of some kind, trading the accuracy of solution for speed. This section presents in short the derivation of this equation and its extension to support light transport in volumes. Next, it explains scattering functions properties. Then, it presents analytical solution to the simplest cases of this equation, and finally, it describes simplifications of light transport equation made by various accelerated rendering algorithms.. 2.2.1. Surface Only Scattering. The following derivation is due to [Veach, 1997]. Global light transport equation is based on a formal definition of local surface scattering. Whenever a beam of light from a direction ω i hits a surface at a point x, it generates irradiance equal to: dE(x, ω i ) = Li (x, ω i )|ω i ◦ N (x)|dσ(ω i ) = Li (x, ω i )dσ ⊥ (ω i ).. (2.7). It can be observed that radiance reflected from a particular point on a surface Lo is proportional to irradiance at that point: dLo (x, ω o ) ∝ dE(x, ω i ). Bidirectional Scattering Distribution Function (BSDF), called scattering function for short, is, by definition, constant of this proportionality: fs (x, ω i , ω o ) =. dLo (x, ω o ) dLo (ω o ) = . dE(x, ω i ) Li (x, ω i )dσ ⊥ (ω i ). Local surface scattering equation is defined as: Z Ls (x, ω o ) = fs (x, ω o , ω i )Li (x, ω i )dσ ⊥ (ω i ).. (2.8). (2.9). Ω. This equation describes how much light is reflected from a surface point x in a particular direction ω o , knowing the incident illumination Li . Total light outgoing from a particular surface point x in a particular direction ω o is a sum of scattered light Ls and emitted light Le : Z Lo (x, ω o ) = Le (x, ω o ) + fs (x, ω o , ω i )Li (x, ω i )dσ ⊥ (ω i ). (2.10) Ω. Incident radiance at a particular surface can be computed using outgoing radiance from another surface: Li (x, ω i ) = Lo (T (x, ω i ), −ω i ). (2.11) The ray casting operator T (x, ω) finds nearest ray-surface intersection point for ray starting from x in direction ω. In order to avoid special case when ray escapes to infinity, the whole scene may be enclosed in a huge, ideally black sphere. The equation 2.11 holds because radiance does not change as light travels along straight lines in vacuum. Substituting 2.11 into 2.10 leads to the final form of rendering equation: Z L(x, ω o ) = Le (x, ω o ) + fs (x, ω o , ω i )L(T (x, ω i ), −ω i )dσ ⊥ (ω i ). (2.12) Ω. The incident radiance Li appears no more in this equation, so the subscript from Lo is dropped. This equation is valid for spectral radiance Lλ as well.. 2.2.2. Volumetric Scattering Extension. Classic rendering equation cannot handle light transport in so-called participating media. These media affect radiance when light travels between surfaces. The assumption of surfaces placed in.

(19) CHAPTER 2. LIGHT TRANSPORT THEORY. 8. vacuum is a good approximation when rays travel on short distances in clean air. It fails, however, in simulation of phenomena such as dusty or foggy atmosphere, emission form fire, or large open environments, where light paths may be many kilometres long. The vacuum version of rendering equation (2.12) can be extended to support volumetric effects by modification of ray casting operator T to account for radiance changes. The extensions are due to [Arvo, 1993], [Pauly, 1999] and [Pharr & Humphreys, 2004]. Participating media may affect ray by increasing or decreasing its radiance while it travels. Increase is due to in-scattering and emission, and decrease due to out-scattering and absorption. The whole participating medium may be described by definition of three coefficients – absorption, emission and scattering in every point of 3D space. Moreover, at every point in which scattering coefficient is larger than zero, there must be also provided a phase function, which performs a role similar to BSDF in classic equation. The emission from a volumetric medium is described by the following expression: dL(x + tω, ω) = Lve (x + tω, ω), dt. (2.13). where Lve (x, ω) is volume emittance. Absorption coefficient is defined as following: dL(x + tω, ω) = −σa L(x + tω, ω). dt. (2.14). Intuitively, the absorption coefficient describes how many times radiance is decreased when light travels a unit distance. The absorption coefficient is measured in [m−1 ], and can take any nonnegative real value. Out-scattering decreases ray radiance in similar way as absorption, but uses scattering coefficient σs instead of σa . The total decrease of ray radiance is expressed by extinction coefficient σe = σa + σs . Fraction of light, which is transmitted between points x and x + sω (beam transmittance), is given by the following formula:  s  Z tr(x, x + sω) = exp − σe (x + tω, ω)dt . (2.15) 0. Beam transmittance has two useful properties: tr(x1 , x2 ). = tr(x2 , x1 ). (2.16). tr(x1 , x3 ). = tr(x1 , x2 )tr(x2 , x3 ).. (2.17). The scattering in the participating medium is described by: Z Lvs (x, ω) = σs (x) fp (x, ω o , ω i )Li (x, ω i )dσ(ω i ).. (2.18). Ω. Radiance added to ray per unit distance due to in-scattering and emission can be expressed as: Lvo (x, ω) = Lve (x, ω) + Lvs (x, ω).. (2.19). Assuming that ray travels through participating medium infinitely long, i.e. it never hits a surface, the total ray radiance change due to participating medium is: Z∞ Li (x, ω) =. tr(x, x + tω)Lvo (x + tω, ω)dt. (2.20). 0. When participating media are mixed with surfaces, similarly as in surface-only rendering, ray casting operator T can be used to found nearest ray-surface intersection. Let y = T (x, ω i ) and.

(20) CHAPTER 2. LIGHT TRANSPORT THEORY. 9. s = kx − yk. The radiance Li incoming at a surface can then be expressed in terms of radiance outgoing from other surface, modified by encountered participating media: Zs Li (x, ω i ) = tr(x, y)Lo (y, −ω i ) +. tr(x, x + tω i )Lvo (x + tω i , ω i )dt. (2.21). 0. Expression 2.21 can be substituted into 2.10, which leads to light transport equation generalized to support participating media mixed with surfaces: Z L(x, ω o ) = Le (x, ω o ) + fs (x, ω o , ω i )· Ω   Zs (2.22) · tr(x, y)L(y, −ω i ) + tr(x, x + tω i )Lvs (x + tω i , ω i )dt dσ ⊥ (ω i ). 0. Similarly as in surface-only version, the subscript from Lo is dropped, and the equation holds for spectral radiance as well. The generalized equation is substantially more complex than the surface-only version, and therefore it may be expected that rendering participating media dramatically hurts performance. Because of that, many existing volumetric rendering algorithms make simplifications of some kind to this general form of volumetric rendering equation.. 2.2.3. Properties of Scattering Functions. The domain of scattering function fs , which describes scattering from all surfaces, is the whole Ω. This function is often defined as a union of simpler functions – reflection fr (Bidirectional Reflection Distribution Function, BRDF) and transmission ft (Bidirectional Transmission Distribution Function, BRDF). Reflection happens when both ω i and ω o are on the same side of a surface. Transmission happens when ω i and ω o are on opposite sides of the surface. Transmission is typically modeled as one two-directional function. Therefore, for reflection only direction of surface normal N is important, while for transmission direction as well as sign of N has to be accounted for. Scattering functions have several important properties. First, to conform the laws of physics, BRDFs must be symmetric, i.e. swapping incident and outgoing directions must not change BRDF value: ∀ω i , ω o ∈ Ω+ fr (ω i , ω o ) = fr (ω o , ω i ). (2.23) However, when surface transmits light, the BTDF typically is not symmetric, but the asymmetry is strictly defined as a function of refraction coefficients of the media on the opposite sides of the surface. The same rule applies to phase functions, but obviously in this case there is no difference between reflection and transmission: ∀ω i , ω o ∈ Ω. fp (ω i , ω o ) = fp (ω o , ω i ).. (2.24). Moreover, all BRDFs and BTDFs (and therefore BSDFs) must be energy conserving, i.e. surfaces cannot reflect more light than they receive: Z ∀ω o ∈ Ω R(ω o ) = fs (ω i , ω o )dσ ⊥ (ω i ) ≤ 1. (2.25) Ω. Analogous relationship holds for phase functions: Z ∀ω o ∈ Ω R(ω o ) = fp (ω i , ω o )dσ(ω i ) = 1.. (2.26). Ω. The equation 2.26 differs from 2.25 in two ways. First, the integration is done with respect to ordinary solid angle, and second, there is strict requirement that phase function is a probability distribution (i.e. it must integrate to one)..

(21) CHAPTER 2. LIGHT TRANSPORT THEORY. 2.2.4. 10. Analytic Solutions. The light transport equation can be solved analytically only in trivial cases, useless in practice. However, these analytical solutions, despite being unable to produce any interesting image, can provide a valuable tool in testing light transport algorithms. Obviously, these tests cannot definitely prove that algorithm which passes them is correct, but nevertheless they provide an aid in removing algorithms’ errors and evaluating their speed of convergence. 1 and constant A very simple test scene is a unit sphere, with constant emission Le (x, ω) ≡ 2π 1 BRDF fr (x, ω i , ω o ) ≡ 2π for each point on the sphere and each direction inside sphere. Since the scene is symmetric (neither geometry nor emission or scattering can break the symmetry), it is easy to verify that radiance L measured along each ray inside the sphere must be identical and equal to 1. This scene can be made a bit more complex with a BRDF which is not necessarily constant, but still with constant reflectivity of R = 0.5. The sphere can be filled with homogenous participating medium with an absorption coefficient σa ≡ 0 and arbitrary scattering coefficient σs . These modifications must not change the result returned by tested light transport algorithms.. 2.2.5. Simplifications. Variety of rendering algorithms, especially those which run in real time, do not attempt to solve full light transport equation. These algorithms may be described theoretically by a (substantially) simplified equations, which are actually solved. In Figure 2.2 there are shown results of simplifications of light transport equation compared with full global illumination solution. The simplified results are given without any additional terms, like ambient light, they just show what simplified equations actually describe. First, one-pass hardware accelerated rasterization implemented in commonly available libraries, such as OpenGL or DirectX, computes only one scattering and supports only point light sources. Therefore rasterization solves the following equation: L(x, ω o ) =. n X. \ fr x, ω o , y\ − yi , i − x Ii x. (2.27). i=1. where n is the number of light sources, yi is the position of ith light source, and Ii is its radiant intensity. Only recently, due to advancements in programmable graphics hardware [Rost & LiceaKane, 2009], the fr function can be arbitrarily complex, and non-point lights can be approximated with reasonable precision. Such dramatic simplifications result in possibility of real time rendering, but at the cost of very poor illumination quality. Classic Whitted ray tracer [Whitted, 1980] (see Section 4.1.2) handles multiple reflections, but only from ideal specular surfaces and for point light sources. This may be seen as replacing integral with a sum of radiances of ideally reflected and transmitted rays: L(x, ω o ) =. n X. \ fr x, ω o , y\ − yi + αL(T (x, ω r ), −ω r ) + βL(T (x, ω t ), −ω t ), i − x Ii x. (2.28). i=1. where ω r is reflected ray direction, ω t is transmitted ray direction, and α and β are color coefficients such that 0 ≤ α < 1, 0 ≤ β < 1 and α + β < 1. An improved approach – Cook’s Distributed Ray Tracing [Cook et al., 1984] computes right hand integral, but only once for area light sources or twice for glossy reflections, so it also does not support global illumination. On the other hand, the radiosity algorithms [Cohen & Wallace, 1993] (see Section 4.1.4) handle multiple reflections, but are limited to matte surfaces only. They solve full light transport equation, but with the assumption that fr ≡ πk , 0 < k < 1. There exist less restrictive radiosity algorithms, however, but they are impractical due to excessive memory consumption. Similarly, the volumetric version of light transport equation is often simplified. In rasterization approach, this equation is typically ignored completely – color of all scene elements exponentially.

(22) CHAPTER 2. LIGHT TRANSPORT THEORY. 11. Figure 2.2: Difference between results of hardware rasterization algorithm (left), classic ray tracing (middle) and full global illumination (right). All guessed lighting terms (e.g. ambient light) were disabled deliberately, to show what given technique actually computes. fade to arbitrarily chosen fog color with the distance from the viewer. This is a very poor approximation, and it cannot simulate majority of volumetric effects, like visible beams of light. However, that is the high price of the ability to run rasterization in real time. On the other hand, physically based rendering systems use less drastic simplifications of volumetric effects. For example, Pharr and Humphreys [Pharr & Humphreys, 2004] implemented single scattering approximation, which seems to be a reasonable trade-off between rendering speed and image quality.. 2.3. Image Formation. This section shows how to create an image from computed radiance. It starts by describing a camera as a device emitting importance. Next, it explains the conversion of light transport equation domain to points only from points and unit direction vectors. Finally, it shows how the equation can be rewritten as integral over space of all light transport paths, which sometimes is easier to solve. These formulae are based on [Veach, 1997], [Pauly, 1999] and [Pharr & Humphreys, 2004], with modification that image is a function defined on real values instead of array of numbers.. 2.3.1. Importance. Cameras may be formally described as devices emitting importance. The importance might be though as hypothetical particles like photons, which propagate from camera, against direction of light flow. Intuitively, tracing importance approximates how important for a rendered image various scene parts are. That is, the importance distribution only is enough to define a camera, and different cameras have different importance distributions. The creation of an image is defined as integral of product of emitted importance (from camera) and radiance (from light sources): Z Z Ii = Wei (xlens , ω)Li (xlens , ω)dσ ⊥ (ω)dA(xlens ), (2.29) Alens. Ω. where the Ii is ith pixel, Wei is its emitted importance and lens is treated as a special surface in the scene – an idealized measuring device – which is able to record radiance while not interfering with light transport. The pixel dependence can be removed in the following way. Let the image domain be a unit square, i.e. [0, 1]2 , and u, v image coordinates: (u, v) ∈ [0, 1]2 . The importance is defined as function of point xlens and direction ω, as well as location on image plane (u, v). The image function I is then evaluated using the expression: Z Z I(u, v) = We (u, v, xlens , ω)Li (xlens , ω)dσ ⊥ (ω)dA(xlens ). (2.30) Alens. Ω.

(23) CHAPTER 2. LIGHT TRANSPORT THEORY. 12. The importance We can be obtained from Wei using a filter function: X We (u, v) = Wei f i (u, v).. (2.31). i. The image equation 2.29 uses importance Wei based on particular filter during image creation process. This seriously limits available image post-processing algorithms. On the other hand, the basic modification in equation 2.30 removes this flaw. The modification seems very simple in theoretical formulation of light transport, but it has significant implications to the described later design of rendering algorithms and their parallelization. Obviously, both these equations are correct for spectral radiance as well. The functional representation of image, however, can cause potential problems when emitted radiance Le described with δ distribution is directly visible. For example, point light sources rendered with pinhole cameras cannot be represented with finite I(u, v) values. Therefore, all algorithms implemented for purpose of this thesis explicitly omit directly visible lights described with δ distributions.. 2.3.2. Integral Formulation. The light transport equation can be reformulated to an integral over all possible light transport paths. The following transformation is based on [Veach, 1997] and [Pharr & Humphreys, 2004]. First, radiance and scattering functions can be expressed in different domains: L(x1 , ω) fs (x2 , ω o , ω i ). = L(x1 → x2 ),. where. = fs (x3 → x2 → x1 ),. ω = x\ 2 − x1 where. ω o = x\ 3 − x2. (2.32) and ω i = x\ 2 − x1 . (2.33). The projected solid angle measure σ ⊥ (ω) is transformed to an area measure A(x) as follows: dσ ⊥ (ω) = V (x1 ↔ x2 ). | cos θ1 || cos θ2 | kx1 − x2 k. 2. dA(x2 ) = G(x1 ↔ x2 )dA(x2 ),. (2.34). where ω = x\ 2 − x1 , θ1 and θ2 are angles between ω and N (x1 ) or N (x2 ), respectively, and V (x1 ↔ x2 ) is the visibility factor – which is equal to 1 if x1 and x2 are mutually visible, and 0 otherwise. Substituting 2.33, 2.32 and 2.34 into 2.12 leads to a light transport equation defined over all scene surfaces instead over all direction vectors: Z L(x2 → x1 ) = Le (x2 → x1 ) + fs (x3 → x2 → x1 )L(x3 → x2 )G(x3 ↔ x2 )dA(x3 ), (2.35) A. and, assuming x0 ≡ xlens and We ≡ We (u, v), image creation equation: Z I(u, v) = We (x0 → x1 )L(x1 → x0 )G(x1 ↔ x0 )dA(x1 )dA(x0 ).. (2.36). A2. Next, using the recursive substitution of the right hand expression of 2.35 into L in 2.36 it is obtained: Z I(u, v) = We (x0 → x1 )G(x1 ↔ x0 )Le (x1 → x0 )dA(x1 )dA(x0 ) + 2 ZA + We (x0 → x1 )G(x1 ↔ x0 )fs (x2 → x1 → x0 ) · A3. · G(x2 ↔ x1 )Le (x2 → x1 )dA(x2 )dA(x1 )dA(x0 ) + Z We (x0 → x1 )G(x1 ↔ x0 )fs (x2 → x1 → x0 )G(x2 ↔ x1 ) ·. + A4. · fs (x3 → x2 → x1 )G(x3 ↔ x2 )Le (x3 → x2 )dA(x3 ) · · · dA(x0 ) + + ··· = ∞ Z X = i=0. Ai+2. We (x0 → x1 )αi Le (xi+1 → xi )dµ(xi ),. (2.37).

(24) CHAPTER 2. LIGHT TRANSPORT THEORY. 13. where αi. i Y. = G(x1 ↔ x0 ). fs (xj+1 → xj → xj−1 )G(xj+1 ↔ xj ). j=1. dµ(xi ). = dA(x0 )dA(x1 ) · · · dA(xi+1 ).. The expressions 2.12, 2.35, and 2.37 for evaluating radiance L are obviously equivalent, but different rendering algorithms often prefer particular form over another. The volumetric version of light transport equation written as integral over light paths was formulated by Pauly [Pauly, 1999]. The main idea of this concept is integration over paths built from any combination of surface and volume scattering. Let bi be the ith bit of binary representation of b, b ∈ N. Let x ¯ = x0 x1 . . . xk be the light transport path with k + 1 vertexes. Integration domain is defined as: A, if bi = 0 k Ψb = (ψ0 ψ1 · · · ψi · · · ψk ), ψi = . (2.38) V, if bi = 1 The integration measure is: µkb (¯ x) = (ψ0 ψ1 · · · ψi · · · ψk ),. ψi =. dA(xi ), if bi = 0 . dV (xi ), if bi = 1. (2.39). The geometric factor is redefined as: Gx (x1 ↔ x2 ) = V (x1 ↔ x2 )tr(x1 ↔ x2 ). c(x1 )c(x2 ). 2,. kx1 − x2 k. c(x) =. | cos θ|, if x ∈ A . 1, if x ∈ V. (2.40). The emitted radiance is Lex (x2 → x1 ) =. Le (x2 → x1 ), if x2 ∈ A , Lev (x2 → x1 ), if x2 ∈ V. (2.41). and scattering function is fx (x2 → x1 → x0 ) =. fs (x2 → x1 → x0 ), if x1 ∈ A . σs fp (x2 → x1 → x0 ), if x1 ∈ V. (2.42). Finally, the integral formulation of volumetric light transport equation is: i+1. I(u, v) =. ∞ 2X −1 Z X i=0. b=1. Ψi+1 2b. We (x0 → x1 )αix Lex (xi+1 → xi )dµi+1 x), 2b (¯. (2.43). where αix is similar to αi , but is defined on Gx and fx instead of G and fs . This equation implicitly makes the common sense assumption that x0 , which is a point on camera lens, is always surface point, i.e. no volumetric sensors are allowed and importance is always emitted from lens surface. However, the radiance emission from volume is accounted properly.. 2.3.3. Image Function. In basic version, image contains only intensity values. However, we found that depth of the first scattering (relative to camera) is potentially very useful in many of postprocessing techniques. Thus, the image function obtained during rendering process is: I : [0, 1]2 × λ −→ R+ ,. (2.44). Id : [0, 1]2 −→ R+ ∪ ∞ ∪ Ve ,. (2.45). where the ∞ means that ray escapes to infinity and Ve means that ray was scattered in participating medium. In the latter case there is no way to reliably define depth of the first intersection..

(25) Chapter 3. Monte Carlo Methods The light transport equation has an analytical solution only for very simple scenes, which are useless in practice, thus appropriate numerical algorithms are necessary to solve it. Classic quadrature rules, e.g. Newton-Cotes or Gaussian quadratures are not well suited to light transport. Light transport integrals are very high-dimensional and integrated functions are discontinuous, which result in poor convergence of quadrature rules based on regular grids. However, usage of algorithm of non-deterministic sampling of integrated functions gives much better results. Non-deterministic algorithms use random numbers to compute the result. According to [Pharr & Humphreys, 2004], they can be grouped in two broad classes – the Las Vegas algorithms, where random numbers are used only to accelerate computations in average case, with final deterministic result (e.g. Quick Sort with randomly selected pivot) and Monte Carlo algorithms, which gives correct results on average. For example, the result of Monte Carlo integration, is not certainly correct, but nevertheless has strict probabilistic bounds on its error. All non deterministic rendering algorithms, from mathematical point of view, can be seen as variants of Monte Carlo integration. The purpose of this chapter is the explanation of statistical methods used in rendering algorithms. Since these methods are well-known, the discussion is brief. Mathematical statistics is explained in detail in [Pluci´ nska & Pluci´ nski, 2000]. The good reference books on general Monte Carlo methods are [Fishman, 1999] and [Gentle, 2003], and [Niederreiter, 1992] on quasiMonte Carlo techniques. Applications of Monte Carlo methods in computer graphics are presented in [Veach, 1997] and [Pharr & Humphreys, 2004]. The rest of this chapter starts with a short review of statistic terms and basic Monte Carlo integration techniques. Next, a distinction between biased and unbiased integration algorithms is explained. This is followed by a description of selected most useful variance (i.e. error) reduction techniques. Finally, some quasi-Monte Carlo methods are presented.. 3.1. Monte Carlo Integration. This section starts with description of concepts employed in statistics. These ideas are then used to construct estimators which approximate integrals of arbitrary functions. Finally, there is a brief analysis of error and convergence rate of Monte Carlo integration, based on variance.. 3.1.1. Statistical Concepts. Let Ψ be a set (called sample space)1 . Let A be a σ-algebra on Ψ and B be the σ-algebra of Borel sets on R. The measure P defined on A is probability if P (Ψ) = 1. The function X : Ψ → R 1 Typically sample space is denoted by Ω. However, in this work Ω is reserved for space of unit directional vectors, thus sample space is denoted by Ψ to avoid confusion.. 14.

(26) CHAPTER 3. MONTE CARLO METHODS. 15. is a single-dimensional random variable if: ∀x ∈ R. X −1 ((−∞, x)) = {ψ : X(ψ) < x} ∈ A.. (3.1). The probability distribution of random variable is defined as: PeX (S) = P ({ψ : X(ψ) ∈ S}) ,. ∀S∈B .. (3.2). The cumulative distribution function (cdf ) is: cdf (x) = PeX ((−∞, x)),. ∀x∈R .. (3.3). Cumulative distribution function cdf (x) may be interpreted as a probability of an event that randomized X value happens to be less or equal given x: cdf (x) = P r {X ≤ x} .. (3.4). The corresponding probability density function (pdf or p) is: pdf (x) =. d cdf (x) . dx. (3.5). Let X1 , X2 , . . . , Xn be random variables and B n be the σ-algebra of Borel sets on Rn . The vector X = (X1 , X2 , . . . , Xn ) is a multidimensional random variable if PeX (S) = P ({ψ : X(ψ) ∈ S}) ,. ∀S∈Bn .. (3.6). The cdf of multidimensional random variable is: cdf (x) = PeX ((−∞, x)),. ∀x∈Rn ,. (3.7). and is corresponding pdf : ∂ n cdf (x) . (3.8) ∂x1 ∂x2 · · · ∂xn The relationship between pdf and cdf can be expressed in more general way using measure theory: Z d cdf (x) and cdf (x) = pdf (x)dµ(x). (3.9) pdf (x) = dµ(x) D pdf (x) =. The expected value of random variable (single- or multidimensional) Y = f (X) is defined as: Z E[Y ] = f (x)pdf (x)dµ(x), (3.10) Ψ. and its variance is:. h i 2 V [Y ] = E (Y − E[Y ]) .. (3.11). Standard deviation σ, which is useful in error estimation, is defined as square root of variance: p (3.12) σ[X] = V [X]. Expected value and variance have the following properties for each α ∈ R: E[αX]. = αE[X],. (3.13). V [αX]. = α2 V [X].. (3.14). The expected value of sum of random variables is a sum of expected values: "N # N X X E Xi = E [Xi ] . i=1. (3.15). i=1. Similar equation holds for variance if and only if the random variables are independent. Using these expressions and some algebraic manipulation variance can be reformulated: h i 2 V [X] = E (X − E[X]) = E X 2 − 2XE[X] + E[X]2 = E X 2 − E[X]2 . (3.16).

(27) CHAPTER 3. MONTE CARLO METHODS. 3.1.2. 16. Estimators of Integrals. Let I be the integral to evaluate: Z I=. f (x)dµ(x).. (3.17). Ψ. The basic Monte Carlo estimator of this integral is: N 1 X f (Xi ) I˜ ≈ FN = , N i=1 pdf (Xi ). (3.18). where ∀x : f (x) 6= 0 pdf (x) > 0. Using definition of expected value it may be found that expected value of estimator 3.18 is equal to integral 3.17: # " Z N N Z 1 X f (x) 1 X f (Xi ) = pdf (x)dµ(x) = E [FN ] = E f (x)dµ(x), (3.19) N i=1 pdf (Xi ) N i=1 Ψ pdf (x) Ψ thus the estimator 3.18 produces correct result on average. Variance of such estimator can be expressed as: # "N # " N N X 1 X 1 1 1 X Xi = 2 V Xi = 2 V [Xi ] = V [F ]. (3.20) V [FN ] = V N i=1 N N N i=1 i=1 This expression for variance is valid if and only if Xi are independent. The variance of single sample V [F ] is equal to: V [F ] = E[F 2 ] − E[F ]2 =. Z Ψ. f 2 (x) dµ(x) − pdf (x). Z. 2 f (x)dµ(x). (3.21). Ψ. Convergence rate of estimator FN can be obtained from Chebyshev’s inequality: ( ) r V [X] P r |X − E[X]| ≥ ≤ δ, δ. (3.22). which holds for any fixed threshold δ > 0 and any random variable X which variance V [X] < ∞. Substitution of estimator FN into Chebyshev’s inequality yields: ( ) r 1 V [F ] P r |FN − I| ≥ √ ≤ δ, (3.23) δ N √ thus for any fixed threshold δ the error decreases with rate O 1/ N. 3.1.3. Biased and Unbiased Methods. All non-deterministic integration algorithms can be divided into two fundamental groups – unbiased and biased. Bias β in statistics is defined as difference between the true value of estimated quantity Q and expected value of the estimator FN : β [FN ] = E [FN ] − Q.. (3.24). Thus, the unbiased algorithms produce correct results on average, without any systematic errors. However, the shortcoming of these algorithms is variance, which is likely to be larger than in biased ones. The variance appears as noise in the rendered image. Non trivial scenes require a lot of computation to reduce this distracting artifact to an acceptable level. On the other hand, biased methods exhibit systematic error. For example, given point in a rendered image can be.

(28) CHAPTER 3. MONTE CARLO METHODS. 17. regularly too bright regardless of number of samples N evaluated. Biased methods still can be consistent (unbiased asymptotically), as long as the error decreases to zero with increasing amount of computation: lim E [FN ] = Q. (3.25) N →∞. The bias is usually difficult to estimate, and even if rendered image does not appear noisy, it still can have substantial inaccuracies, typically seen as excessive blurring or illumination artifacts on fine geometrical details. However, despite non-zero bias, biased methods tends to converge substantially faster (they have lower variance) than unbiased ones.. 3.2. Variance Reduction Techniques. The basic Monte Carlo estimator potentially can have high variance, which directly leads to poor efficiency. In this section there is a brief revision of general methods which substantially reduce variance of this estimator, without introducing bias. On the other hand, biased methods used in computer graphics are often dedicated to particular light transport algorithms, and therefore are described in Chapter 4.. 3.2.1. Importance and Multiple Importance Sampling. A variance of Monte Carlo estimator 3.18 can be decreased, when the pdf (x) is made more proportional to f (x). Intuitively, this approach tends to place relatively more samples wherever integrand is large, therefore reducing integration error. Particularly, when pdf (x) ∝ f (x) is satisfied exactly, the variance is zero: −1. Z pdf (x) = cf (x),. c=. f (x)dµ(x). ,. (3.26). Ψ. and therefore Z V [F ] = Ψ. f 2 (x) dµ(x) − pdf (x). Z. 2 Z 2 Z f (x)dµ(x) = c−1 f (x)dµ(x) − f (x)dµ(x) = 0.. Ψ. Ψ. Ψ. However, in order to achieve zero variance, the function must be integrated analytically, to obtain normalization constant c. This is impossible, since otherwise there would not be the necessity of using Monte Carlo integration at all. Fortunately, using a pdf which is proportional, or almost proportional, to at least one of factors of f , typically decreases variance. This technique is called Importance Sampling. Suppose that f (x) can be decomposed into product f (x) = f1 (x)f2 (x) · · · fn (x) and there exist probability densities pdfi proportional (or roughly proportional) to each factor fi . If the standard Importance Sampling is used, the pdfi used for sampling f (x) has to be chosen at algorithm design time. This can have disastrous consequences for algorithm efficiency, if the chosen pdfi poorly matches the overall f (x) shape. In this case, Importance Sampling can actually increase variance over sampling with uniform probability. In this case Multiple Importance Sampling [Veach & Guibas, 1995] can be used. This technique was designed to improve the reliability of Importance Sampling when the appropriate pdf cannot be chosen at the design time. The main idea of this method is to define more than one pdf (each of them is potentially good candidate for importance sampling) and let the algorithm chose the best one at runtime, when the actual shape of integrand is known. The algorithm does this by computing appropriate weights and returning the estimator as weighted sum of samples from these pdf s: Fnm =. m n X f (Xij ) 1 X wi (Xij ) , m j=1 pdfi (Xij ) i=1. ∀x. n X i=1. wi (x) = 1.. (3.27).

(29) CHAPTER 3. MONTE CARLO METHODS. 18. The appropriate choice of weights wi is crucial for obtaining low variance estimator. According to [Veach & Guibas, 1995] following set of weights is the optimal choice: pdfi (x) . wi (x) = Pn j=1 pdfj (x). (3.28). However, Multiple Importance Sampling causes an issue that has large impact on design and implementation of sampling routines. The standard Importance Sampling technique requires just two methods – sampling points xi with given probability pdf (xi ) and evaluating f (xi ). The Multiple Importance Sampling, however, requires an additional operation – computing pdf (x) for an arbitrary argument x. Intuitively, the operation may be interpreted as ‘compute the hypothetical probability of returning given x’. This operation is usually more difficult than computing probability during sampling xi , since algorithm has no knowledge of random choices necessary to select arbitrary value, as it has in sampling procedure. Fortunately, when approximated probabilities are computed, the Multiple Importance Sampling estimator still is correct, but crude approximation hurts performance.. 3.2.2. Russian Roulette and Splitting. Russian roulette and Splitting are designed to adaptively change sampling density without introducing bias. These two techniques were introduced to computer graphics by Arvo and Kirk [Arvo & Kirk, 1990]. Suppose that estimator F is sum of estimators F = F1 + F2 + . . . + Fn . Russian roulette allows random skipping evaluation of these terms: ( Fi −(1−qi )c , with probability qi , 0 qi (3.29) Fi = c, with probability 1 − qi , where c is arbitrary constant, typically c = 0. When the estimator is a sum of infinite number of terms S = F1 + F2 + . . . + Fi + . . ., the Russian roulette still can be applied, provided that the sum S is finite. Let Si = Fi + Fi+1 + . . . be the partial sum. Then S can be reexpressed as S = S1 , S1 = F1 + S2 , S2 = F2 + S3 , . . .. Each sum Si is then evaluated with probability qi , and set to 0 otherwise. Provided that at most finite number of qi is equal to 1 and ∃ ε > 0 : qi < 1 − ε, for almost all qi , evaluation of sum S is randomly terminated with probability 1 after n terms. This leads to expression:    X i n Y 1 1 1 1   Fi  S0 = (3.30) F1 + F2 + . . . + Fn . . . = q1 q2 qn q j i=1 j=1 Russian roulette, however, increases variance of the estimator. Nevertheless, since it reduces computational cost of it, Russian roulette can improve the estimator efficiency (product of variance and cost), if probabilities qi are chosen carefully. Moreover, Russian roulette can be used to terminate computation of infinite series without introducing statistical bias. Splitting works in opposite direction to Russian roulette. Splitting increases the number of samples in order to reduce variance. Splitting increases computation time, but nevertheless, if performed carefully, can improve sampling efficiency. According to [Veach, 1997], the splitting technique works as follows: n 1X Fij , (3.31) Fi0 = n i=1 where Fij are independent samples from Fi .. 3.2.3. Uniform Sample Placement. Random samples tend to clump together leaving large portions of domain relatively empty. Uneven distribution of samples leads to increased variance of estimators and therefore low sampling.