0032137: Coding Rules - merge redundant .lxx files into header files within Package gp
[occt.git] / dox / specification / pbr_math.md
CommitLineData
6b6d06fa 1PBR math (rasterization) {#specification__pbr_math}
4e8371cb 2========================
3@tableofcontents
4
27ce2519 5@section pbr_preface Preface
4e8371cb 6
7**Empirical** illumination models like **Phong reflection model** have been used in real-time graphics for a long time due to their simplicity, convincing look and affordable performance.
8Before programmable pipeline has been introduced, graphics cards implemented Gouraud shading as part of fixed-function Transformation & Lighting (T&L) hardware blocks.
9Nowadays, however, numerous trade-offs of this simplicity (like lighting partially baked into object material properties and others) pushed developers to **Physically-Based Rendering** (**PBR**) illumination models.
10
11PBR models try to fit surface shading formulas into constrains of physical laws of light propagation / absorption / reflection - hence, called "physically-based".
12There are two main categories of PBR illumination:
13
14 1. Non-real-time renderer (cinematic).
15 2. Real-time renderer.
16
17The main objective of cinematic renderer is uncompromised quality, so that it relies on ray-tracing (path-tracing) rendering pipeline.
18Although performance of current graphics hardware does not make it possible using computationally-intensive path-tracing renderer in real-time graphics, it can be used in interactive fashion.
19
20"Physically-based" does not necessarily mean physically-correct/precise.
21The main objective of real-time PBR renderer is to be fast enough even on low-end graphics hardware.
22So that in contrast, it hardly relies on rasterization rendering pipeline, various approximations and tricks making it applicable in real-time, while looking good enough and preserving some physical properties.
23
24OCCT 3D Viewer provides both kinds of PBR renderers, and although they share some details in common, this article is devoted to real-time PBR metallic-roughness illumination model.
25This article describes the math underneath PBR shading in OCCT 3D Viewer and its GLSL programs.
26However, this article does not clarifies related high-level APIs nor PBR material creation pipelines, as this is another topic.
27
27ce2519 28@section pbr_notation Notation
4e8371cb 29
30| | | |
31|-:|:-|:-|
32| \f$n\f$ | normal (on surface) | \f$\|n\|=1\f$ |
33| \f$v\f$ | view direction | \f$\|v\|=1\f$ |
34| \f$l\f$ | light | \f$\|l\| = 1\f$ |
35| \f$h=\frac{v+l}{\|v + l\|}\f$ | half vector | |
36| \f$m\f$ | metallic factor | \f$[0, 1]\f$ |
37| \f$r\f$ | roughness factor | \f$[0, 1]\f$ |
38| \f$IOR\f$ | index of refraction | \f$[1, 3]\f$ |
39| \f$c\f$ | albedo color | \f$(R, G, B)\f$ |
40
41\f$\cos\theta_l=(n \cdot l)\f$
42
43\f$\cos\theta_v=(n \cdot v)\f$
44
45\f$\cos\theta_h=(n \cdot h)\f$
46
47\f$\cos\theta_{vh}=(v \cdot h)\f$
48
27ce2519 49@section pbr_illumination_model Illumination model
4e8371cb 50
51The main goal of illumination model is to calculate outgoing light radiance \f$L_o\f$ along the certain direction.
52The starting point of calculation might be the view direction \f$v\f$ aimed from point on surface (or in more general case just in space) to viewer position.
53Considering the point on opaque surface with normal \f$n\f$ the main equation of illumination can be defined as:
54
55\f[L_o=\int\limits_H f(v, l) L_i(l) \cos\theta_l\, \mathrm{d}l\f]
56
57Where \f$L_i(l)\f$ is light radiance coming from \f$l\f$ direction, \f$f(v,l)\f$ is **Bidirectional Reflectance Distribution Function** (**BRDF**) and \f$H\f$ is hemisphere which is oriented regarding to the surface normal \f$n\f$.
58Opaqueness of the surface mentioned earlier is important because in that case hemisphere is enough.
59More general model will require to consider directions all around a whole sphere and is not observed in this paper.
60\f$\cos\theta_l\f$ factor appearing is caused by affection of surface area and light direction mutual orientation to the amount of radiance coming to this area.
61This is mainly due to geometric laws. The rest part of integral is the key of the whole illumination model.
62BRDF defines it's complexity and optical properties of material.
a84f14ab 63It has to model all light and material interactions and also has to satisfy some following criteria in order to be physical correct [@ref Duvenhage13]:
4e8371cb 64* Positivity: \f$f(v,l) \geq 0\f$
65* Helmholtz reciprocity: \f$f(v,l) = f(l, v)\f$ (follows from 2<sup>nd</sup> Law of Thermodynamics)
66* Energy conservation: \f$\displaystyle \forall v \, \int\limits_H f(v,l) \cos\theta_l \, \mathrm{d}l = 1\f$ (in order not to reflect more light than came)
67
68It is worth to be mentioned that \f$f(v,l)\f$ depends on \f$n\f$ also but it is omitted to simplify notation. BRDF is usually split into two parts:
69
70\f[f(v,l) = f_d(v,l)+f_s(v, l)\f]
71
72Where \f$f_s(v, l)\f$ (specular BRDF) models reflection light interaction on surface and \f$f_d(v,l)\f$ (diffuse BRDF) models other processes happening depth in material (subsurface scattering for example).
73So that illumination equation might be rewritten as:
74
75\f[L_o=\int\limits_H (f_d(v,l)+f_s(v, l)) L_i(l) \cos\theta_l\, \mathrm{d}l\f]
76
a84f14ab 77PBR theory is based on **Cook-Torrance specular BRDF** [@ref Cook81]. It imagines surface as set of perfectly reflected micro faces distributed on area in different ways which is pretty good model approximation of real world materials.
4e8371cb 78If this area is small enough not to be able to recognize separate micro surfaces the results becomes a sort of averaging or mixing of every micro plane illumination contribution.
79In that level it allows to work with micro faces in statistical manner manipulating only probabilities distributions of micro surfaces parameters such as normals, height, pattern, orientation etc.
80In computer graphics pixels are units of images and it usually covers a relatively large areas of surfaces so that micro planes can be considered to be unrecognizable.
81Going back to the BRDF the Cook-Torrance approach has the following expression:
82
83\f[f_s(v,l)=\frac{DGF}{4\cos\theta_l\cos\theta_v}\f]
84
85Three parts presented in nominator have its own meaning but can have different implementation with various levels of complexity and physical accuracy.
86In that paper only one certain implementation is used. The \f$D\f$ component is responsible for **micro faces normals distribution**.
87It is the main instrument that controls reflection's shape and strength according to **roughness** \f$r\f$ parameter.
a84f14ab 88The implementation with good visual results is **Trowbridge-Reitz GGX** approach used in Disney's RenderMan and Unreal Engine [@ref Karis13]:
4e8371cb 89
90\f[D=\frac{\alpha^2}{\pi(\cos^2\theta_h(\alpha^2-1) + 1)^2}\f]
91
92Where \f$\alpha = r^2\f$. This square in needed only for smoother roughness parameter control.
93Without it the visual appearance of surface becomes rough too quickly during the parameter's increasing.
94
95The second \f$G\f$ component is called **geometric shadowing** or attenuation factor.
a84f14ab 96The point is that micro surfaces form kind of terrain and can cast shadows over each other especially on extreme viewing angles [@ref Heitz14].
97**Smith-Schlick model** [@ref Heitz14], [@ref Schlick94] has been chosen as implementation:
4e8371cb 98
99\f[\displaystyle G=\frac{\cos\theta_l \cos\theta_v}{(\cos\theta_l(1-k)+k)(\cos\theta_v(1-k)+k)}\f]
100
101Where \f$k=\frac{\alpha}{2}\f$, which means \f$k=\frac{r^2}{2}\f$ in terms of this paper.
102But \f$G\f$ depends on many factors so that it's approximations has float nature and can be modified a little bit in some cases in order to get more pleasant visual results.
103One of this modification will be described later in following chapters.
104
105The last component \f$F\f$ shows **how much light is reflected from surface** and is called **Fresnel's factor**.
106The rest amount of radiance might be absorbed or refracted by material.
107The most accurate expression of it is pretty complicate for calculation so that there is a variety of approximations.
a84f14ab 108The good one with less computation efforts is **Schlick's implementation** [@ref Schlick94]:
4e8371cb 109
110\f[F=F_0+(1-F_0)(1-\cos\theta_{vh})^5\f]
111
112Here \f$F_0\f$ is material's response coefficient at normal incidence (zero angle).
113Fresnel's factor has to be calculated differently for metals and dielectric/non-metals, but PBR theory tries to come up with universal formula for all types of material.
114In order to do that it is needed to be noticed that Schlick's approximation is applicable only to non-conductors and in that case \f$F_0 = F_{dielectric} = \left(\frac{1-IOR}{1+IOR}\right)^2\f$.
115**Index of Refraction** \f$IOR\f$ shows the proportion between light speed in vacuum (or even in air) and in material.
116The reference value of \f$IOR\f$ for plastic is **1.5**, and this value can be considered as default for all unknown dielectrics.
117In practice this parameter controls reflectance ability of material.
a84f14ab 118Also it should be remembered that this approximation produces poor results with large \f$IOR\f$ values so that it is recommended to be kept in range of \f$[1, 3]\f$ in order to get plausible Fresnel's factor [@ref Lazanyi05], [@ref Lagarde13].
4e8371cb 119This formula might be further propagated onto metals by using \f$F_0\f$ measured specifically for certain metal.
120It can be considered as some kind of a 'color' of metal and can be stored as albedo parameter \f$c\f$.
121And the final step of defining Fresnel's factor formula is mixing all this \f$F_0\f$ using metallic parameter \f$m\f$ (**metalness**):
122
123\f[F_0 = F_{dielectric}(1-m)+cm\f]
124
125For pure dielectrics with \f$m=0\f$ exactly Schlick's approximation will be used.
126For pure metals with \f$m=1\f$ it will be a little inaccurate but the same formula with measured \f$F_0\f$ values.
127Everything else for \f$m \in (0, 1)\f$ is not physically correct and it is recommended to keep \f$m\f$ exactly 1 or 0.
128Intermediate values may represent mixed areas for smooth transition between materials - like partially rusted metal (rust is mostly dielectric).
129Also it might be useful when parameters are read from textures with filtering and smoothing.
130
131BRDF described above has one important trait making computations easier called **isotropy**.
132Isotropy in this case means independence from rotation about normal resulting from supposition of uniform micro faces distribution at any direction along a surface.
133It allows to simplify random samples generation during Monte-Carlo integrals calculation and reduce dimensions of some lookup tables, which will be discussed in following chapters.
134Of course, isotropic materials form only subset of all real world's materials, but this subset covers majority of cases.
135There are special models considering special anisotropic traits of surfaces like a grinding of metal or other with dependency on rotation about normal;
136these models require special calculation tricks and additional parameters and are out of scope of this paper.
137
138The only thing left to do is to define \f$f_d(v,l)\f$.
139This part is responsible for processes happening in depth of material.
140First of all the amount of input light radiance participating in these processes is needed to be calculated.
141And it exactly can be realized from already known Fresnel's factor \f$F\f$ showing amount of reflected light but in negative term in this case in order to get the radiance left after reflection:
142
143\f[1-F\f]
144
145This part of ingoing light is assumed to be refracted in depth of surface and variety of events may happen there.
146A sequence of absorptions, reflections and reemissions more or less leads to light's subsurface scattering.
147Some part of this scattered light can go back outside but in modified form and in pretty unpredictable directions and positions.
148For opaque materials this part is noticeable and forms it's own color.
149If subsurface's paths of light are small enough and points of output are distributed locally around the input point it's possible to work in statistical way similar to the micro faces.
150This assumption covers a big amount of real world opaque materials.
151Other materials like skin, milk etc. with noticeable effect of subsurface scattering usually presented in form of partial translucency and some kind of self emission
152have more widely distributed output points and require more accurate and complicate ways of modeling with maybe some theory and techniques from volumetric rendering.
153The simple but visually enough assuming for statistically driven type of materials is just the same radiance for any direction. It results to **Lambertian's BRDF**:
154
155\f[\frac{c}{\pi}\f]
156
157Where \f$\pi\f$ is normalization coefficient in order to meet BRDF's criteria and \f$c\f$ is material's own color formed by adventures of light under surface.
158There is one detail about light interaction bringing some physicality to the model, and that is an absence of this diffuse component in metals.
159Metals reflect main part of light and the rest of it is absorbed being transformed into other form (mostly heat).
160That is the main visual difference between metallic and non-metallic materials realizing of which brings model to higher level of quality in compare to older non-physical models.
161
162So that all parts described above can be combined into united diffuse BRDF:
163
164\f[f_d(v,l) = (1-F)(1-m)\frac{c}{\pi}\f]
165
166\f$m\f$ is recommended to be exactly 1 or 0 but all values between can represent transition areas, as mentioned before.
167
168In this chapter one possible implementation of illumination model reflecting main PBR principles has been defined.
169The next step is using of it in practice.
170
27ce2519 171@section pbr_practical_application Practical application
4e8371cb 172
173It's time to apply deduced illumination model in practice.
174And the first step of it is separation of **direction based light sources** from illumination integral.
175Directional nature of such light sources means possibility to calculate it's influence to point of surface using only one direction and its intensity.
176Usually sources of this type do not have physical size and are represented only by position in space (for point or spot lights) or by direction itself (direction light imagined to be too far point sources like sun).
177This is just a kind of abstraction, while real world light emitters have noticeably sizes.
178But sources with realistic form and size cannot be presented in discrete term and require continuous integrals calculations or special approximations in order to be accurately injected to the model.
179In most cases direct based light sources in form of emitting points in space or just certain directions are good approximations and are enough for beginning.
180Having finite discrete amount of it in scene and considering only single direction from every of these lights, the integral is transformed just to the sum:
181
182\f[L_{direct} = \sum_{j=1}^M f(v, l_j) L_i^{direct}(l_j) \cos\theta_{l_j}\f]
183
184Where \f$M\f$ is a number of sources, \f$l_j\f$ is a direction and \f$L_i^{direct}\f$ is an intensity related to this direction.
185\f$direct\f$ label means that illumination has been computed directly from sources.
186The BRDF can be used directly without any calculation problems.
a84f14ab 187The only exception might be \f$k\f$ in \f$G\f$ factor - it is recommended to use \f$ k = (r+1)^2 / 8 \f$ in order to get more pleasant results [@ref Karis13] (that is modification mentioned in previous chapter).
4e8371cb 188And actually it is enough to finally see something.
189There will be correct visualization with assumption of complete dark environment and absence of other points influence.
190It is called **local illumination**. Based on this name there is also a global or **indirect illumination** and that is the rest of integral:
191
192\f[L_{indirect} = \int\limits_H f(v, l) L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l\f]
193
194It includes influence of light reflected or scattered from other points and environment's contribution.
7b5e532f 195It's impossible to achieve photorealistic results without this component, but it is also very difficult to compute.
4e8371cb 196While the cross point light interaction cannot be calculated in a simple way (especially in real time rendering), the environment illumination has some options to be realized via precomputational work before visualization.
197But right now lets summarize the practical application of illumination model.
198At this moment the output radiance is represented as:
199
200\f[L_o = L_{direct} + L_{indirect}\f]
201
202Where \f$L_{direct}\f$ is direction based light sources contribution which can be directly computed just applying bare formulas.
203It is enough to get some results in terms of local illumination but without \f$L_{indirect}\f$ component image will not be looked lifelike.
204\f$L_{indirect}\f$ is not trivial thing for calculation and that is stumbling block for real time rendering applications.
205But it can be relatively easy implemented in case of environment illumination via some precomputational work about which will be told in details in following chapters.
206
27ce2519 207@section pbr_image_based_lighting Image based lighting
4e8371cb 208
209The next goal after \f$L_{direct}\f$ calculation is to find \f$L_{indirect}\f$.
210And it would be easier if \f$L_i^{indirect}(l)\f$ was known for every \f$l\f$.
211That is the main assumption of **image based lightning** (**IBL**).
212In practice, it can be achieved using environment image map, which is a special image representing illumination from all possible directions.
213This image might be a photo capturing a real world environment (spherical 360 degrees panoramas) or generated image baking the 3D scene itself, including in that case reflections of other objects.
214Environment image might be packed in different ways - **cube maps** and equirectangular maps are the most commonly used.
215Anyway, it allows \f$L_i^{indirect}(l)\f$ to be defined for every \f$l\f$ and its practical implementation in form of images gives name for this approach.
216Lets back to indirect illumination integral:
217
218\f[L_{indirect} = \int\limits_H f(v, l) L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l\f]
219
220Substituting the BRDF by its expression allows to split indirect illumination into diffuse and specular components:
221
222\f[L_{indirect} = \int\limits_H f_d(v,l)L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l + \int\limits_H f_s(v,l)L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l = \f]
223
224\f[= (1-m)\frac{c}{\pi}\int\limits_H (1-F)L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l + \frac{1}{4}\int\limits_H \frac{DFG}{\cos\theta_l \cos\theta_v}L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l\f]
225
226This splitting seems not to lead to simplicity of calculation but these two parts will be computed in slightly different ways in future.
227Lets write down this separately:
228
229\f[L_{indirect}^d = (1-m)\frac{c}{\pi}\int\limits_H (1-F)L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l\f]
230
231\f[L_{indirect}^s = \frac{1}{4}\int\limits_H \frac{DFG}{\cos\theta_v \cos\theta_l} L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l\f]
232
233Next transformations of these expressions require understanding of numerical way to find hemisphere integral and also its performance optimization techniques.
234And that the topic of the next chapter.
235
27ce2519 236@section pbr_monte_carlo_integration Monte-Carlo numeric integration
4e8371cb 237
238**Monte-Carlo** is one of numeric methods to **find integral**.
239It is based on idea of mathematical expectation calculation.
240In one dimensional case if \f$f(x)\f$ is a function with parameter distributed according to probability density \f$p(x)\f$ the mathematical expectation of it can be found using following expression:
241
242\f[E = \int\limits_{-\infty}^\infty f(x) p(x)\, \mathrm{d}x\f]
243
244Also this expectation can be approximated in statistical term using certain sequence of random variable \f$x\f$:
245
246\f[E \approx \frac{1}{N} \sum_{i=1}^{N} f(x_i)\f]
247
248It can be used in general definite integrals calculations.
249Just valid \f$p(x)\f$ defined on \f$[a, b]\f$ range and sequence \f$x_i\f$ generated according to it are needed for that:
250
251\f[\int\limits_a^b f(x)\, \mathrm{d}x = \int\limits_a^b \frac{f(x)}{p(x)}p(x)\, \mathrm{d}x = \int\limits_{-\infty}^{\infty} \frac{f(x)}{p(x)}p(x)\, \mathrm{d}x \approx \frac{1}{N} \sum_{i=1}^{N} \frac{f(x_i)}{p(x_i)}\f]
252
253Where \f$f(x)\f$ is considered to be equal to zero outside of \f$[a, b]\f$ range.
254This is also true for functions on sphere or hemisphere:
255
256\f[\int\limits_{H|S} f(l)\, \mathrm{d}l \approx \frac{1}{N}\sum_{i=1}^{N} \frac{f(l_i)}{p(l_i)}\f]
257
258The main questions are choosing \f$p(l)\f$ and generating samples \f$l_i\f$.
259The one of the simple ways is uniform distribution along sphere or hemisphere.
260Lets realize that on sphere for example.
261There are \f$4\pi\f$ possible directions in terms of sphere's areas and steradians (direction can be presented as dot on a unit sphere):
262
263\f[\int\limits_S 1\, \mathrm{d}l = 4\pi\f]
264
265Where \f$S\f$ is the unit sphere.
266In order to be uniform \f$p(l)\f$ must be constant and satisfy normalization criteria:
267
268\f[\int\limits_S p(l)\, \mathrm{d}l = 1\f]
269
270So that \f$p(l) = \frac{1}{4\pi}\f$.
271Usually direction \f$l\f$ is parametrized by spherical coordinates \f$\phi \in [0, 2\pi]\f$ and \f$\theta \in [0, \pi]\f$ boiling down to the 2D samples generation.
272But in these terms joint \f$p(\theta, \phi)\f$ will be looked slightly different due to variables transition.
273\f$l\f$ is defined in regular Cartesian coordinates \f$l=(x, y, z)\f$ with \f$\|l\| = 1\f$.
274The spherical coordinates transform looks like:
275
276\f[x = r\sin\theta\cos\phi,\, y = r\sin\theta\sin\phi,\, z = r\cos\theta\f]
277
278Where \f$r = 1\f$.
279In order to express probability density using new variables it is needed to multiply this density by Jacobian of transform:
280
281\f[p(\theta,\phi) = p(l)|J_T|\f]
282
283Where:
284
285\f[|J_T| = \begin{vmatrix} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \phi} \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} & \frac{\partial y}{\partial \phi} \\ \frac{\partial z}{\partial r} & \frac{\partial z}{\partial \theta} & \frac{\partial z}{\partial \phi} \end{vmatrix} = r^2\sin\theta\f]
286
287So that joint probability density in new variables looks like:
288
289\f[p(\theta, \phi) = \frac{\sin\theta}{4\pi}\f]
290
291This variable transfer rule of **Probability Density Function** (**PDF**) will be useful in following chapters, when integral calculation optimization techniques will be being told.
292Having \f$p(\theta, \phi)\f$ the partial single dimensional probability densities are able to be found:
293
294\f[p(\phi) = \int\limits_0^\pi p(\theta, \phi)\, \mathrm{d}\theta = \frac{1}{4\pi} \int\limits_0^\pi \sin\theta\, \mathrm{d}\theta = \frac{1}{2\pi}\f]
295
296\f[p(\theta) = \int\limits_0^{2\pi} p(\theta, \phi)\, \mathrm{d}\phi = \frac{\sin\theta}{4\pi}\int\limits_0^{2\pi}1\, \mathrm{d}\phi = \frac{\sin\theta}{2}\f]
297
298The final step is sequence generation itself.
299In order to be able to generate values with arbitrary distributions it is helpful to start from uniform numbers in range of \f$[0, 1]\f$.
300And that can be done via any known true- and pseudo- random generators.
301Even simple \f$\frac{1}{i}\f$ sequence is appropriate for beginning but it can be not so efficient in terms of computations convergence.
302There are specially designed series for the last reason and it will be tackled in chapter about optimizations.
303The \f$\phi\f$ variable is noticed to be uniformly distributed so that it can be directly generated without any additional manipulations.
304Just range \f$[0, 1]\f$ is needed to be mapped to range \f$[0, 2\pi]\f$.
305For any other variables including \f$\theta\f$ the inverse transform sampling approach can be applied.
306First of all **cumulative distribution function** (**CDF**) is needed to be found.
307It is probability of random value to be less than argument of this functions by definition.
308For continues distributions it can be expressed in following form:
309
310\f[F(x) = \int\limits_{-\infty}^x p(x')\, \mathrm{d}x'\f]
311
312Lets find CDF for \f$\theta\f$:
313
314\f[F(\theta) = \int\limits_{-\infty}^\theta p(\theta')\, \mathrm{d}\theta' = \int\limits_0^\theta \frac{\sin\theta'}{2}\, \mathrm{d}\theta' = \frac{1-\cos\theta}{2}\f]
315
316The CDF maps \f$\theta\f$ values from range of \f$[0, \pi]\f$ to probability in range of \f$[0, 1]\f$.
317The next step is to find inverse cumulative function which can be not so trivial sometimes but pretty obvious in current case:
318
319\f[F^{-1}(u) = \arccos(1-2u)\f]
320
321If substitute uniform distributed in range \f$[0, 1]\f$ values \f$u\f$ as argument of this function the values with origin probability density will appear.
322In other words:
323
324\f[\theta = \arccos(1 - 2u),\, u \in [0, 1],\, p(u) = 1 \Rightarrow p(\theta) = \frac{\sin\theta}{2}\f]
325
326That is the key of this random values generation technique.
327All steps described above can be also done for hemisphere:
328
329\f[p(l) = \frac{1}{2\pi}\f]
330
331\f[p(\theta, \phi) = \frac{\sin\theta}{2\pi}\f]
332
333\f[p(\phi) = \int\limits_0^\frac{\pi}{2} p(\theta, \phi)\, \mathrm{d}\theta = \frac{1}{2\pi} \int\limits_0^\frac{\pi}{2} \sin\theta\, \mathrm{d}\theta = \frac{1}{2\pi}\f]
334
335\f[p(\theta) = \int\limits_0^{2\pi} p(\theta, \phi)\, \mathrm{d}\phi = \frac{\sin\theta}{2\pi}\int\limits_0^{2\pi}1\, \mathrm{d}\phi = \sin\theta\f]
336
337\f[\theta = \arccos(1-u)\f]
338
339Mote-Carlo integration cannot guarantee exact estimation of convergence speed with using random generated samples.
340There is only probability estimation of it.
341But this algorithm is pretty universal and relatively simple to be applied to almost any function using at least uniform distributed points.
342Moreover special \f$p(l)\f$ can be chosen and special pseudo-random sequences can be designed in order to speed up convergence for some functions (following chapter talk about that in details).
343That is why this method is widely used in computer graphics and demonstrates good results.
344Also another one advantage is worth to be mentioned - possibility to iteratively accumulate computations and present intermediate results during rendering which is used in some ray tracing applications.
345
27ce2519 346@section pbr_split_sum Split sum
4e8371cb 347
348Lets go back to the image based lighting and the figure of specular component.
349As was defined before that is hemisphere integral with following expression:
350
351\f[L_{indirect}^s = \frac{1}{4}\int\limits_H \frac{DFG}{\cos\theta_v \cos\theta_l} L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l\f]
352
353The Monte-Carlo integration algorithm can be directly applied:
354
355\f[L_{indirect}^s = \int\limits_H f_s(v, l)L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l \approx \frac{1}{N}\sum_{i=1}^N \frac{f_s(v, l_i) L_i^{indirect}(l_i) \cos\theta_{l_i}}{p(v, l_i)}\f]
356
357\f$p(v, l_i)\f$ depends on \f$v\f$ and implicitly on \f$r\f$ in order to be completely general.
358Optimization strategies use different samples distributions for different view direction orientations and roughness values.
359Anyway even with all optimization techniques this algorithm continues to require too much calculations.
360Good visual results require noticeable number of samples and using this approach for every point in real time rendering becomes unrealistic.
361The way to avoid these enormous calculations is doing them beforehand somehow.
a84f14ab 362The first trick on the way to this is split the sum separating environment light component [@ref Karis13]:
4e8371cb 363
364\f[L_{indirect}^s \approx \frac{1}{N} \sum_{i=1}^N \frac{f_s(v, l_i) L_i^{indirect}(l_i) \cos\theta_{l_i}}{p(v, l_i)} \approx \left( \frac{1}{N} \sum_{i=1}^N L_i^{indirect}(l_i) \right) \left( \frac{1}{N} \sum_{i=1}^N \frac{f_s(v, l_i) \cos\theta_{l_i}}{p(v, l_i)} \right)\f]
365
366Where the second brackets represent approximation of integral so that the expression can be rewritten as:
367
368\f[L_{indirect}^s \approx \frac{1}{N} \sum_{i=1}^N \frac{f_s(v, l_i) L_i^{indirect}(l_i) \cos\theta_{l_i}}{p(v, l_i)} \approx \left( \frac{1}{N} \sum_{i=1}^N L_i^{indirect}(l_i) \right) \int\limits_H f_s(v, l) \cos\theta_l\, \mathrm{d}l\f]
369
370This integral is exact \f$L_{indirect}^s\f$ in condition when \f$L_i^{indirect}(l) = 1\f$ what just means white uniform environment.
371The sum before it is kind of averaged environment illumination.
372The main accomplishment after all this manipulations is possibility to calculate light and BRDF components separately.
373The sum with \f$L_i^{indirect}(l_i)\f$ can be computed beforehand for every normal direction and stored to image called specular map but with some valuable details.
374The problem is that \f$l_i\f$ samples must be generated according to \f$p(v, l_i)\f$ distribution depended on \f$v\f$ and \f$r\f$ as was mentioned earlier.
375Variation of normal is not enough in that case and these variables are needed to be considered too.
376The ways to resolve it are topic of one of the following chapters and now understanding the fact that at least this part can be precomputed before rendering is enough for now.
377And it is important not to miss out that there is no more BRDF influence in this sum and only \f$p(v, l)\f$ can affect in this case.
378That is why it is so important to strict to PDF during samples generation and that is why \f$p(v, l)\f$ must be correlated with BRDF somehow in this approximation approach with splitting.
379For example completely mirroring materials with \f$r = 0\f$ will not be looked as expected if just uniform distribution is used
380because such surfaces have only one possible direction from which light can be reflected along view direction in compare with \f$N\f$ absolutely scattered in case of uniform or many other distributions.
381
382The rest part also can be saved to image. Lets unroll its expression:
383
384\f[\int\limits_H f_s(v, l) \cos\theta_l\, \mathrm{d}l = \int\limits_H \frac{DGF}{4\cos\theta_v \cos\theta_l} \cos\theta_l\, \mathrm{d}l\f]
385
386This integral is not actually a scalar.
387That is RGB value due to only \f$F\f$ factor and even more only to \f$F_0\f$.
388In order to simplify future computations \f$F_0\f$ is needed to be moved out of integral.
389Substitution of Schlick's approximation helps to achieve it:
390
391\f[F = F_0+(1-F_0)(1-\cos\theta_{vh})^5 = F_0(1-(1-\cos\theta_{vh})^5) + (1-\cos\theta_{vh})^5\f]
392
393\f[\int\limits_H \frac{DGF}{\cos\theta_v \cos\theta_l} \cos\theta_l\, \mathrm{d}l = F_0 \int\limits_H \frac{DG}{4\cos\theta_v \cos\theta_l} (1-(1-\cos\theta_{vh})^5) \cos\theta_l\, \mathrm{d}l + \int\limits_H \frac{DG}{4\cos\theta_v \cos\theta_l} (1-\cos\theta_{vh})^5 \cos\theta_l\, \mathrm{d}l\f]
394
395This form may not look easier, but it has several advantages.
396The first one is independence from globally defined \f$L_i^{indirect}(l)\f$, so that normal orientation does not matter and can be set in any handful way for calculations (Z axis for example).
397The second one results from isotropic illumination system allowing \f$\phi\f$ component of view vector to be set arbitrarily (0 for example) and \f$\cos\theta_v\f$ will be enough to define view direction.
398And the third one is scalar nature of integrals so that only two precomputed numbers are needed to find BRDF part of \f$L_{indirect}^s\f$.
399Considering dependency of these integrals from \f$\cos\theta_v\f$ and \f$r\f$ both of it can be precomputed and stored to 2D look-up image variating these two parameters in range \f$[0, 1]\f$ with two channels consisting of scale and bias for \f$F_0\f$.
400
401Current result for \f$L_{indirect}^s\f$ is computing it using 2D image for BRDF part and omnidirectional image for environment illumination.
402There were a lot of words about Monte-Carlo optimizations techniques and about PDF choice which is important not only in terms of numeric integration but in terms of visual results correctness.
403It's time to talk about that.
404
27ce2519 405@section pbr_importance_sampling Importance sampling
4e8371cb 406
407Current goal is to speed up Monte-Carlo integration of Cook-Torrance like integrals with following expression:
408
409\f[\int\limits_H \frac{DG}{4\cos\theta_v \cos\theta_l} g(v, l) \cos\theta_l\, \mathrm{d}l\f]
410
411Where \f$g(v, l)\f$ is just arbitrary function representing Fresnel's factor itself or its components.
412In order to increase convergence the samples with larger contribution (or in other words with larger function's values) have to appear more frequently than others proportionally to its contribution.
413So that less significant summand with less influence to result will be considered rarely and in opposite way parts brining noticeable changes to the sum will be taken often.
414That is the main idea of **importance sampling technique**.
415\f$p(l)\f$ has to represent significance of sample in terms of integrated function via probability somehow.
416And it will be like that if PDF is already part of original function because in that case probability density directly affects to contribution forming.
417Separating this distribution component is one possible and effective way to realize importance sampling strategy.
418In integral presented above PDF part already exists and that is \f$D\f$ component.
419But it is distribution of micro faces normals or ideally reflection direction or \f$h\f$ in other word and not light directions distribution which is needed in fact.
420Anyway that is good starting point and lets generate \f$h\f$ vectors first.
421\f$D\f$ has the following expression:
422
423\f[D=\frac{\alpha^2}{\pi(\cos^2\theta_h(\alpha^2-1) + 1)^2}\f]
424
425Frankly speaking \f$D(h)\f$ is called normal distribution but cannot be directly used as hemisphere distribution.
426Originally it is statistical factor used to define total area of micro faces \f$\mathrm{d}A_h\f$
a25d5aaa 427whose normals lie within given infinitesimal solid angle \f$\mathrm{d}h\f$ centered on \f$h\f$ direction using the original small enough area of macro surface \f$\mathrm{d}A\f$ [@ref Walter07]:
4e8371cb 428
429\f[dA_h = D(h)\,\mathrm{d}h\, \mathrm{d}A\f]
430
431First of all this factor must be positive:
432
433\f[D(h) \geq 0\f]
434
435But the total area of micro faces landscape is at least equal to origin surface but even bigger in general:
436
437\f[1 \leq \int\limits_H D(h)\, \mathrm{d}h\f]
438
439This trait does not allow to use \f$D\f$ as hemisphere distribution.
440But it can be fixed with following feature:
441
442\f[\forall v\, \int\limits_H D(h)(v \cdot h)\, \mathrm{d}h = (v \cdot n)\f]
443
444Which means that total area of micro faces projected to any direction must be the same as projected area of origin macro surface.
445It is pretty tricky trait in \f$D\f$ definition but it leads to interesting results in condition when \f$v = n\f$:
446
447\f[\int\limits_H D(h)\cos\theta_h\, \mathrm{d}h = 1\f]
448
449So that \f$\cos\theta_h\f$ coefficient normalizes normal distribution in terms of hemisphere and allows to use it as distribution.
450Finally PDF of half vectors can be wrote:
451
452\f[p(\theta_h, \phi_h) = D\cos\theta_h\sin\theta_h = \frac{\alpha^2 \cos\theta_h\sin\theta_h}{\pi(\cos^2\theta_h(\alpha^2-1) + 1)^2}\f]
453
454\f$\sin\theta_h\f$ results from spherical coordinate system transfer which was described in Monte-Carlo integration chapter.
455Lets strict to samples generation procedure and find partial probability densities:
456
457\f[p(\phi_h) = \int\limits_0^\frac{\pi}{2} p(\theta_h, \phi_h)\, \mathrm{d}\theta_h = \int\limits_0^\frac{\pi}{2} \frac{\alpha^2 \cos\theta_h\sin\theta_h}{\pi(\cos^2\theta_h(\alpha^2-1) + 1)^2}\, \mathrm{d}\theta = \frac{1}{2\pi}\f]
458
459\f[p(\theta_h) = \int\limits_0^{2\pi} p(\theta_h, \phi_h)\, \mathrm{d}\phi_h = \int\limits_0^{2\pi} \frac{\alpha^2 \cos\theta_h\sin\theta_h}{\pi(\cos^2\theta_h(\alpha^2-1) + 1)^2}\, \mathrm{d}\phi = \frac{2 \alpha^2 \cos\theta_h\sin\theta_h}{(\cos^2\theta_h(\alpha^2-1) + 1)^2}\f]
460
461\f$p(\phi_h)\f$ is unnecessary to be calculated analytically.
462The fact of independency from \f$\phi\f$ is enough to figure out that this coordinate is uniformly distributed.
a84f14ab 463Anyway the \f$F(\theta_h)\f$ is next step [@ref Cao15]:
4e8371cb 464
465\f[F(\theta_h) = \int\limits_0^{\theta_h} \frac{2 \alpha^2 \cos\theta'_h\sin\theta'_h}{(\cos^2\theta'_h(\alpha^2-1) + 1)^2}\, \mathrm{d}\theta'_h = \int\limits_{\theta_h}^0 \frac{2 \alpha^2}{(\cos^2\theta'_h(\alpha^2-1) + 1)^2}\, \mathrm{d}(\cos^2\theta'_h) = \frac{\alpha^2}{\alpha^2-1}\int\limits_0^{\theta_h} \mathrm{d}\frac{1}{\cos^2\theta'_h(\alpha^2-1)+1} =\f]
466
467\f[ = \frac{\alpha^2}{\alpha^2-1} \left( \frac{1}{\cos^2\theta_h(\alpha^2-1) + 1} - \frac{1}{\alpha^2} \right) = \frac{\alpha^2}{\cos^2\theta_h(\alpha^2-1)^2+(\alpha^2-1)} - \frac{1}{\alpha^2-1}\f]
468
469In order to apply inverse transform sampling the \f$F^{-1}(u)\f$ is needed to be found:
470
471\f[F^{-1}(u) = \theta_h = \arccos\sqrt{\frac{1-u}{u(\alpha^2-1)+1}}\f]
472
473So that there is no more obstacles to generate \f$h\f$.
474But the main goal was \f$l\f$ direction.
475In order to get it the view vector \f$v\f$ has to be reflected related to already found \f$h\f$:
476
477\f[l = 2(v \cdot h)h - v\f]
478
479That is practical side of light direction generation.
480But the theoretical one is needed to be resolved to calculate sum.
481It is time to find \f$p(l)\f$ using known \f$p(h)\f$.
482First of all the fact that \f$l\f$ is just transformed \f$h\f$ is needed to be understood.
483In that way the light direction's PDF has following expression:
484
485\f[p(l) = p(h)|J_T|\f]
486
487Where \f$|J_T|\f$ is Jacobian of reflection transformation.
488Lets find it.
489Right now \f$n\f$ is axis from which \f$\theta\f$ spherical coordinate is encountered.
490The first step is setting \f$v\f$ as starting point of \f$\theta\f$ instead of \f$n\f$.
491This is linear transform so that \f$|J_T| = 1\f$.
492Next step is transfer to spherical coordinate with \f$|J_T| = \sin\theta_{vh}\f$.
493Due to previous step \f$\theta_{vh}\f$ is used instead of \f$\theta_h\f$.
494In this coordinate system reflecting of \f$v\f$ relative to \f$h\f$ is just doubling \f$\theta_{vh}\f$ and Jacobian of it is equal to \f$\frac{1}{2}\f$.
495In series of transform the Jacobians are multiplied so that currently \f$|J_T| = \frac{1}{2}\sin\theta_{vh}\f$.
496And the final step is inverse transform to Cartesian coordinate system with \f$|J_T| = (\sin\theta_{vl})^{-1} = (\sin2\theta_{vh})^{-1}\f$.
a84f14ab 497Combining this all together the following expression is obtained for reflection transform Jacobian [@ref Schutte18]:
4e8371cb 498
499\f[|J_T| = \frac{\sin\theta_{vh}}{2\sin2\theta_{vh}} = \frac{\sin\theta_{vh}}{4\sin\theta_{vh}\cos\theta_{vh}} = \frac{1}{4\cos\theta_{vh}}\f]
500
501And finally \f$p(l)\f$ looks like:
502
503\f[p(l) = p(h)|J_T| = \frac{D\cos\theta_h}{4\cos\theta_{vh}}\f]
504
505Lets go back to the Monte-Carlo sum and insert found result to it:
506
507\f[\int\limits_H \frac{DG}{4\cos\theta_v \cos\theta_l} g(v, l) \cos\theta_l\, \mathrm{d}l \approx \frac{1}{N} \sum_{i=1}^N \frac{DG\, g(v, l_i) \cos\theta_{l_i}}{4\cos\theta_v \cos\theta_{l_i}\, p(l_i)} = \frac{1}{N} \sum_{i=1}^N \frac{G\, g(v, l_i) \cos\theta_{l_i} \cos\theta_{vh_i}}{\cos\theta_v \cos\theta_{l_i} \cos\theta_{h_i}}\f]
508
509Here \f$G\f$ component is recommended to be calculated with original \f$k=\frac{\alpha}{2} = \frac{r^2}{2}\f$ in order to get more plausible result.
510Of course, all \f$\cos\f$ must be clamped to range \f$[0, 1]\f$ because integral is calculated on a hemisphere and all expressions are defined on it.
511\f$\cos\theta_v \cos\theta_{l_i}\f$ in denominator can be reduced with exactly the same part in geometric attenuation factor in order to avoid additional zero division cases.
512
513Summarizing importance sampling strategy described above the convergence of Monte-Carlo integration can be improved using special PDF correlated with integrated function.
514In case of BRDF with normal distribution functions \f$D\f$ the PDF producing procedure is defined.
515Practically half vector \f$h\f$ is generated first and \f$l\f$ is obtained from it by view vector \f$v\f$ reflecting.
516Due to this transformation final form of probability density used in sum is quite different but also has defined algorithm of calculation.
517For isotropic Cook-Torrance BRDF the \f$\cos\theta_v\f$ and roughness \f$r\f$ are enough to start generation so that all integrals of that kind can be precalculated in 2D look-up tables variating these two parameters.
518The same samples generation procedure must be used in specular map baking described in next chapter.
519
27ce2519 520@section pbr_specular_map Specular map
4e8371cb 521
522The situation with BRDF part of \f$L_{indirect}^s\f$ is clear now and \f$L_i^{indirect}(l)\f$ sum is left to be discussed.
523That was called **specular map** and has following form:
524
525\f[\frac{1}{N}\sum_{i=1}^N L_i^{indirect}(l_i)\f]
526
527As was mentioned this sum must be calculated for every normal direction using the same samples generation principles as in numeric integration computation.
528This principles require two scalar parameters \f$\cos\theta_v\f$ and \f$r\f$ but now \f$\phi\f$ really matters.
529So that in fact the specular map has to be saved in 3D table consisting omnidirectional textures.
530That is a big expense of computational and memory resources.
531A couple of tricks helps to reduce dimensions.
532First of all the \f$\cos\theta_v\f$ and \f$\phi\f$ can be just excluded.
533In that way \f$v\f$ is considered to be equal to \f$n\f$.
534Of course this approach produces an error and affects the final result.
a84f14ab 535It can be fixed more or less by \f$\cos\theta_{l_i}\f$ weighting [@ref Karis13]:
4e8371cb 536
537\f[\frac{1}{N} \sum_{i=1}^N L_i^{indirect}(l_i) \cos\theta_{l_i}\f]
538
539It is not a complete solution but practice shows that it is enough to get plausible illumination with sacrificing of lengthy reflections at grazing angles which exist in fact if everything is honestly computed.
540The problem is that for \f$v \neq n\f$ considering this sum to be defined related to \f$n\f$ became incorrect.
541For example, for complete mirroring materials with \f$r = 0\f$ this sum must boil down to \f$L_i^{indirect}(v_r)\f$
542but not to \f$L_i^{indirect}(n)\f$ where \f$v_r\f$ is just reflected \f$v\f$ or in other words \f$v_r = 2(v \cdot n)n - v\f$.
543That it just mirroring reflection principle.
544Assumption of \f$n = v\f$ also means that \f$n = v = v_r\f$.
545In that way radiance map better to be considered as averaging of illumination coming from \f$v_r\f$.
546So that it has become to be defined related to reflection direction which has to be calculated before map's fetching.
547
548Anyway, there are just two dimensions in radiance look-up table remain.
549The rest one with \f$r\f$ parameter cannot be reduced.
550There is no other ways except just roughness variation but in order to simplify that computations can be done for several values and the rest ones lying between can be obtained from linear interpolation.
551This is another source of visual artifacts but it also works good in practice and that is pretty common approach.
552But it still requires noticeably amount of samples and that is for every pixel related to each \f$r\f$ value.
553It can be appropriate for precomputations but still limits using dynamic environments in real time rendering or just even static environments but on weak devices such as mobile ones.
554And there are several possible ways to improve radiance map baking performance.
555
556The first one is using textures with smaller resolutions for larger roughnesses.
557The point is that smaller \f$r\f$ values produce map saving more details from origin environment in opposite to larger ones representing lower frequency components and working as low pass filters in fact.
558So less pixels in combination with linear interpolation is enough to store less detailed convolutions.
559Moreover, this approach naturally works with textures levels of details in graphics API
560so that every certain radiance map related to certain \f$r\f$ can be stored on its own mip level and be directly fetched with linear interpolation not only over one texture but over levels too.
561As practice shows 6 levels are enough.
562
563After reducing pixels count it is turn for samples number.
564And again correlation with roughness can be noticed.
565For example map for completely mirroring materials with \f$r = 0\f$ the same sample \f$l_i = v_r\f$ will be produced.
566So that only one sample is enough in this case.
567In opposite way directions for \f$r = 1\f$ will be scattered over almost whole hemisphere what requires as much samples as available.
568The 'locality' of distribution is decreased during increasing roughness and it is possible to assume that samples number might to be proportional to this 'locality' keeping accuracy at the same level.
569But how can 'locality' be interpreted in terms of probability distribution? One possible way is CDF meaning.
570\f$F(\theta_h)\f$ has been already defined and by definition it shows the probability of random value \f$\theta_h\f$ to be less than argument of CDF.
571In other words \f$F(\theta'_h) = u\f$ means that probability of \f$\theta_h\f$ to be in range of \f$[0, \theta'_h]\f$ is \f$u\f$.
572The inverse task of range searching using given \f$u\f$ can be solved with help of \f$F^{-1}(u) = \theta'_h\f$.
573If \f$u\f$ is close to 1 (exact 1 has no sense because in that case \f$\theta'_h = \max\theta_h = \frac{\pi}{2}\f$)
574then \f$\theta'_h\f$ represents the range of the most probable or most frequently generated values and that can be interpreted as 'locality' of distribution.
575After that if samples number of the worst case with \f$r = 1\f$ is set (\f$N(1) = \max N\f$) the other ones can be estimated using following formula:
576
577\f[N(r) = N(1)\frac{\theta'_h(r)}{\frac{\pi}{2}} = N(1)\frac{2\theta'_h(r)}{\pi} = N(1)\frac{2F^{-1}(u)}{\pi} = N(1)\frac{2}{\pi}\arccos\sqrt{\frac{1-u}{u(\alpha^2-1)+1}}\f]
578
579It is approximate expression representing only estimated general proportionality so that cases of \f$r = 0\f$ and \f$r = 1\f$ must be processed separately with \f$N(0) = 1\f$ and \f$N(1) = \max N\f$.
580\f$u\f$ can be parameter of this optimization strategy controlling speed of samples reducing in order to balance performance and quality (\f$u = 1\f$ disables this optimization at all).
581This pretty tricky technique allows reducing calculations for every pixels without sacrificing the quality.
582
583In addition to optimizations mentioned before another one can be applied in order to help to reduce samples numbers as previous one.
584Using less samples produces image noise due to discrete nature of Monte-Carlo approximation.
585But it can be slightly smoothed using some prefiltration.
586The idea is that for the directions with small PDF or in other words for rare directions the samples near of it is unlikely to be generated.
587So that this direction represents the averaged illumination from relatively big area on hemisphere but approximate it by just a constant.
588It wold be better to get from such direction already averaged over bigger area environment.
589It can be achieved using mip levels of origin \f$L_i^{indirect}\f$ whose pixels of one level is exact 4 averaged pixels from previous one.
590Also mip levels generation is build in most common graphic API so there are no problems with it.
591But first of all the area covered by one sample is needed to be found.
a84f14ab 592And that can be done as [@ref Colbert07]:
4e8371cb 593
594\f[\Omega_s = \frac{1}{N\,p(l)} = \frac{4\cos\theta_{vh}}{ND\cos\theta_h}\f]
595
596Circumstance of \f$v = v_r = n\f$ leads to \f$\cos\theta_{vh}\f$ and \f$\cos\theta_h\f$ reducing so expression becomes even simpler:
597
598\f[\Omega_s =\frac{4}{ND}\f]
599
600Of course all zero divisions must be avoided by clamping, for example.
601After that the area covered by one pixel of environment map is calculated.
602In case of a cube map it looks like:
603
604\f[\Omega_p = \frac{4\pi}{6k^2}\f]
605
606Where \f$k\f$ is size of one cube map side in pixels (sides are assumed to be quads).
607Finally the mip level of origin environment map which is needed to be fetched for this certain sample is defined by following expression:
608
609\f[lod = \frac{1}{2} \log_2\left(\frac{\Omega_s}{\Omega_p}\right)\f]
610
611The mathematics connected with mip levels sizes lie behind it but this is out of scope of this paper.
612In combination with previous optimization technique this approach allows \f$N(1)\f$ to be smaller keeping visual results good.
613
614That is not all possible optimization tricks but at least these three significantly reduces compute efforts and brings radiance map calculation to weak devices or even to dynamic environments in real time but in reasonable limits.
615
616In that way \f$L_{indirect}^s\f$ can be completely computed without any integral approximations.
617Only 2D look-up table of BRDF part and mip mapped omnidirectional texture of irradiance map are needed.
618The first one can be got even without any environment.
619It was achieved using some rough approximations and assumptions but despite of that the visual result are still plausible and can be compared even with ray traced images.
620In order to complete whole image based lighting the \f$L_{indirect}^d\f$ component is left to be discussed.
621
27ce2519 622@section pbr_spherical_harmonics Spherical harmonics
4e8371cb 623
624Lets go back to diffuse indirect illumination component represented by following formula:
625
626\f[L_{indirect}^d = (1-m)\frac{c}{\pi}\int\limits_H (1-F)L_i^{indirect}(l)\cos\theta_l\, \mathrm{d}l\f]
627
628Of course, Monte-Carlo algorithm can be applied directly and hemisphere integral can be precalculated for every normal direction
629but dependence from \f$v\f$ in Fresnel's factor does not allow to do it efficiently (every \f$v\f$ direction is needed to be considered again).
27ce2519 630In order to resolve it modified version of Schlick's approximation has been created [[?](TODO)]:
4e8371cb 631
632\f[F \approx F_{ss}=F_0+(\max(1-r, F_0))(1-\cos\theta_v)^5\f]
633
634It differs from origin one and loses accuracy a little bit but now there is no light direction inside
635so that it can be considered as kind of screen space defined Fresnel's factor (\f$ss\f$ means exactly 'screen space') and can be removed from integral:
636
637\f[L_{indirect}^d = (1-m)(1-F_{ss})\frac{c}{\pi} \int\limits_H L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l\f]
638
639The resulted expression without \f$(1-m)\f$ and \f$(1-F_{ss})\f$ parts is pretty known entity called **irradiance**.
640It can be precalculated using \f$\cos\theta_l\f$ as PDF for importance sampling (actually it is only option in this case excluding uniform distribution).
641But even with that samples will be scattered almost over whole hemisphere.
642As was discussed in previous chapter this case requires significant amount of samples in order to average illumination with appropriate quality.
643Poor accuracy resulted from lack of summand can be noticed especially on high frequency environments having a lot of contrasting details.
644It worth to be mentioned that irradiance is used not only in BRDF.
645Omnidirectional diffuse illumination captured for certain point or even for several points uniformly or hierarchically distributed is base of some baking global illumination techniques.
646There it is called a **light probe**. So that other way to compute and store irradiance maps was found resolving many mentioned problems.
647The Fourier's decomposition analogue for spherical function allows to achieve this.
648That would be easy to explain concept directly on example.
649So lets start from \f$L_i^{indirect}(l)\f$.
650It is spherical function because directions are just points on sphere.
651The decomposition looks like:
652
653\f[L_i^{indirect}(l) = \sum_{i = 0}^\infty \sum_{j=-i}^i L_i^j y_i^j(l)\f]
654
655Where \f$y_i^j(l)\f$ are spherical functions forming orthonormalized basis called spherical harmonics and \f$L_i^j\f$ is decompositions coefficients.
656Orthonormality means that dot product of two basis elements is equal to 1 if this is the same functions and is equal to zero otherwise.
657Dot product on a sphere is defined as integral of functions multiplication. In other words:
658
659\f[\int\limits_S y_i^j(l)\, y_{i'}^{j'}(l)\, \mathrm{d}l = \begin{cases} 1 & \quad i,j = i',j' \\ 0 & \quad \mathrm{otherwise}\end{cases}\f]
660
a84f14ab 661Function basis with such traits is known and is described by following formulas [@ref Guy18]:
4e8371cb 662
663\f[y_i^{j > 0}(\theta, \phi) = \sqrt{2}K_i^j\cos(j\phi)P_i^j(\cos\theta)\f]
664\f[y_i^{j<0}(\theta, \phi) = \sqrt{2}K_i^j\sin(j\phi)P_i^{|j|}(\cos\theta)\f]
665\f[y_i^0(\theta, \phi) = K_i^0P_i^0(\cos\theta)\f]
666
667\f[K_i^j = \sqrt{\frac{(2i+1)(i-|j|)!}{4\pi(i+|j|)!}}\f]
668\f[P_0^0(x) = 1\f]
669\f[P_1^0(x) = x\f]
670\f[P_i^i(x) = (-1)^i(2i-1)!!(1-x^2)^\frac{i}{2}\f]
671\f[P_i^j(x) = \frac{(2i-1)xP_{i-1}^j(x) - (i + j - 1)P_{i-2}^j}{i - j}\f]
672
673Here \f$K_i^j\f$ are normalization factors and \f$P_i^j\f$ are **Legendre's polynomials**.
674Decomposition coefficients \f$L_i^j\f$ are dot product of origin function (\f$L_i^{indirect}(l)\f$ in current case) and corresponding basis element. It can be written down as:
675
676\f[L_i^j = \int\limits_S L_i^{indirect}(l)\,y_i^j(l)\, \mathrm{d}l\f]
677
678Fact that all calculation happen over a sphere but not over hemisphere right now is important not to be missed.
679That was example of spherical function decomposition but not a solution for original task which looks like:
680
681\f[\int\limits_H L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l\f]
682
683First of all, lets transform this integral to be defined not over hemisphere but sphere:
684
685\f[\int\limits_H L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l = \int\limits_S L_i^{indirect}(l)\overline{\cos}\theta_l\, \mathrm{d}l\f]
686
687Where \f$\overline{\cos}\f$ is cosine clamped to zero which can be expressed as:
688
689\f[\overline{\cos}\theta_l = \max(\cos\theta_l, 0)\f]
690
691Resulted expression can be considered as convolution in terms of spherical functions where \f$L_i^{indirect}(l)\f$ is target and \f$\overline{\cos}\theta_l\f$ is core.
692This integral may seem independent but in fact hemisphere is oriented related to \f$n\f$ therefore \f$\overline{\cos}\theta_l\f$ depends on it too and became a kind of 'oriented' version of cosine.
693That is pretty tricky and explanation about meaning of convolution on sphere is out of scope of this paper.
a84f14ab 694Fact that this is convolution analogue related to \f$n\f$ is enough for now [@ref Aguilar17], [@ref Ramamoorthi01].
4e8371cb 695The goal of looking at integral from this angle is using of convolution's trait allowing to compute decomposition using just only coefficients of function and core.
696\f$\overline{\cos}\theta_l\f$ is independent from \f$\phi_l\f$ and in case of such radial symmetric cores the resulting coefficients boil down to following formula:
697
698\f[(L_i^{indirect}(l) \ast \overline{\cos}\theta_l)_i^j = \frac{1}{K_i^0}L_i^j\, c_i^0 = \sqrt{\frac{4\pi}{2i+1}}L_i^j\, c_i^0\f]
699
700Where \f$c_i^0\f$ are spherical harmonics factors corresponding to \f$\overline{\cos}\theta\f$.
701Therefore the final decomposition looks like:
702
703\f[\int\limits_{H(n)} L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l = \int\limits_S L_i^{indirect}(l)\overline{\cos}\theta_l\, \mathrm{d}l = \sum_{i=0}^\infty \sum_{j = -i}^i \sqrt{\frac{4\pi}{2i+1}}L_i^j\, c_i^0\, y_i^j(n)\f]
704
705\f$c_i^0\f$ is left to be found.
706Due to independence from \f$\phi\f$ all \f$c_i^{j \neq 0} = 0\f$.
707The rest ones are calculated by regular dot product with basis functions:
708
709\f[c_i^0 = c_i = \int\limits_S y_i^0(l)\, \overline{\cos}\theta_l\, \mathrm{d}l = \int\limits_0^{2\pi} \mathrm{d}\phi \int\limits_0^\pi y_i^0(\theta, \phi)\, \overline{\cos}\theta \sin\theta\, \mathrm{d}\theta = \int\limits_0^{2\pi} \mathrm{d}\phi \int\limits_0^\frac{\pi}{2} y_i^0(\theta, \phi)\, \cos\theta\sin\theta\, \mathrm{d}\theta = \f]
710
711\f[= 2\pi\int\limits_0^\frac{\pi}{2} y_i^0(\theta, \phi)\, \cos\theta\sin\theta\, \mathrm{d}\theta = 2\pi K_i^0\int\limits_0^\frac{\pi}{2} P_i^0(\cos\theta)\, \cos\theta\sin\theta\, \mathrm{d}\theta\f]
712
713\f$\sin\theta\f$ appears due to transfer from integral over sphere to double integral where \f$\mathrm{d}l = \sin\theta\, \mathrm{d}\theta\, \mathrm{d}\phi\f$.
714There is an analytical solution for this expressions:
715
716\f[c_1 = \sqrt{\frac{\pi}{3}}\f]
717
718\f[c_{odd} = 0 \quad c_{even} = 2\pi\sqrt{\frac{2i+1}{4\pi}}\frac{(-1)^{\frac{i}{2}-1}}{(i+2)(i-1)}\frac{i!}{2^i\left(\frac{i!}{2}\right)^2}\f]
719
720Starting from about the third \f$c_i\f$ the coefficients become less and less valuable so that only couple of them is enough in order to approximate \f$\overline{\cos}\theta\f$ with appropriate accuracy.
721The same principle is true for convolution too because its coefficients are multiplied by \f$c_i\f$.
722So there is no need to use more than even three bands (\f$i = 0, 1, 2\f$) of basis functions.
a84f14ab 723Lets write down them all in Cartesian coordinate [@ref Ramamoorthi01]:
4e8371cb 724
725\f[y_0^0 = \frac{1}{2}\sqrt{\frac{1}{\pi}} = Y_0^0\f]
726
727\f[y_1^{-1} = -\frac{1}{2}\sqrt{\frac{3}{\pi}}y = Y_1^{-1}y\f]
728\f[y_1^0 = \frac{1}{2}\sqrt{\frac{3}{\pi}}z = Y_1^0z\f]
729\f[y_1^1 = -\frac{1}{2}\sqrt{\frac{3}{\pi}}x = Y_1^1x\f]
730
731\f[y_2^{-2} = \frac{1}{2}\sqrt{\frac{15}{\pi}}xy = Y_2^{-2}xy\f]
732\f[y_2^{-1} = -\frac{1}{2}\sqrt{\frac{15}{\pi}}yz = Y_2^{-1}yz\f]
733\f[y_2^0 = \frac{1}{4}\sqrt{\frac{5}{\pi}}(3z^2-1) = Y_2^0(3z^2-1)\f]
734\f[y_2^1 = -\frac{1}{2}\sqrt{\frac{15}{\pi}}xz = Y_2^1xz\f]
735\f[y_2^2 = \frac{1}{4}\sqrt{\frac{15}{\pi}}(x^2-y^2) = Y_2^2(x^2-y^2)\f]
736
737All \f$Y_i^j\f$ are just constants so that it can be moved from integral during calculations and can be taken from precomputed table.
738Other constants related to \f$c_i\f$ can be united and also be calculated beforehand:
739
740\f[\hat{c}_i = \frac{1}{K_i^0}\, c_i = \sqrt{\frac{4\pi}{2i+1}}\, c_i\f]
741
742Finally expression of irradiance map approximation can be defined:
743
744\f[\int\limits_{H(n)} L_i^{indirect}(l) \cos\theta_l\, \mathrm{d}l \approx \sum_{i=0}^2 \sum_{j=-i}^i L_i^j\, \hat{c}_i\, y_i^j(n)\f]
745
746Where \f$\hat{c}_i\f$ is precalculated constants \f$y_i^j(n)\f$ are pretty easy functions and only \f$L_i^j\f$ are needed to be precomputed.
747Of course \f$L_i^j\f$ are integrals over even whole sphere but now there is only nine of it instead of one for every pixel of omnidirectional image.
748Moreover, texture is not needed at all in that case and only 9 colors representing \f$L_i^j\f$ have to be saved.
749The Monte-Carlo algorithm can be applied with just uniform samples distribution without importance sampling at all.
750\f$Y_i^j\f$ are used twice: in \f$L_i^j\f$ calculations and in sum after that.
751So there is sense to store only squares of it.
a84f14ab 752All tables with constants presented below [@ref Ramamoorthi01]:
4e8371cb 753
754| |
755|-|
756| \f$(Y_0^0)^2 \approx (0.282095)^2\f$ |
757| \f$(Y_1^{-1})^2 = (Y_1^0)^2 = (Y_1^1)^2 \approx (0.488603)^2\f$ |
758| \f$(Y_2^{-2})^2 = (Y_2^{-1})^2 = (Y_2^1)^2 \approx (1.092548)^2\f$ |
759| \f$(Y_2^0)^2 \approx (0.315392)^2\f$ |
760| \f$(Y_2^2)^2 \approx (0.546274)^2\f$ |
761
762| | |
763|-|-|
764| \f$\hat{c}_0\f$ | \f$3.141593\f$ |
765| \f$\hat{c}_1\f$ | \f$2.094395\f$ |
766| \f$\hat{c}_2\f$ | \f$0.785398\f$ |
767
768Summarizing all mathematics above spherical harmonics decomposition boils down irradiance map to only 9 values which is needed to be precalculated as integrals.
769As practice shows this is very good approximation of diffuse indirect illumination component.
770
27ce2519 771@section pbr_transparency Transparent materials
4e8371cb 772
773TODO
774
27ce2519 775@section pbr_low_discrepancy Low discrepancy sequence
4e8371cb 776
777TODO
27ce2519 778
779@section pbr_references References
780
a84f14ab 781<table cellpadding="4">
782<tr><td valign="top">
783@anchor Duvenhage13 **[Duvenhage13]**
784</td><td>
27ce2519 785Bernardt Duvenhage, Kadi Bouatouch, D.G. Kourie,
786"Numerical Verification of Bidirectional Reflectance Distribution Functions for Physical Plausibility",
787*Proceedings of the South African Institute for Computer Scientists and Information Technologists Conference*,
788October 2013.
a84f14ab 789</td></tr>
27ce2519 790
a84f14ab 791<tr><td valign="top">
792@anchor Cook81 **[Cook81]**
793</td><td>
27ce2519 794Robert Cook, Kenneth Torrance,
795"A Refectance Model for Computer Graphics",
796*SIGGRAPH '81: Proceedings of the 8th annual conference on Computer graphics and interactive techniques*,
797August 1981, pp. 307-316.
a84f14ab 798</td></tr>
27ce2519 799
a84f14ab 800<tr><td valign="top">
801@anchor Karis13 **[Karis13]**
802</td><td>
27ce2519 803Brian Karis, "Real Shading in Unreal Engine 4", *SIGGRAPH 2013 Presentation Notes*.
a84f14ab 804</td></tr>
27ce2519 805
a84f14ab 806<tr><td valign="top">
807@anchor Heitz14 **[Heitz14]**
808</td><td>
27ce2519 809Eric Heitz, "Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs",
810*Journal of Computer Graphics Techniques*, Vol. 3, No. 2, 2014.
a84f14ab 811</td></tr>
27ce2519 812
a84f14ab 813<tr><td valign="top">
814@anchor Schlick94 **[Schlick94]**
815</td><td>
27ce2519 816Christophe Schlick, "An inexpensive brdf model for physically-based rendering",
817*Computer Graphics Forum 13*, 1994, pp. 233-246.
a84f14ab 818</td></tr>
27ce2519 819
a84f14ab 820<tr><td valign="top">
821@anchor Lazanyi05 **[Lazanyi05]**
822</td><td>
27ce2519 823Istvan Lazanyi, Lazslo Szirmay-Kalos, "Fresnel term approximations for Metals", January 2005.
a84f14ab 824</td></tr>
27ce2519 825
a84f14ab 826<tr><td valign="top">
827@anchor Lagarde13 **[Lagarde13]**
828</td><td>
27ce2519 829Sebastien Lagarde, "Memo on Fresnel equations",
830*Blog post*: [https://seblagarde.wordpress.com/2013/04/29/memo-on-fresnel-equations/](https://seblagarde.wordpress.com/2013/04/29/memo-on-fresnel-equations/).
a84f14ab 831</td></tr>
27ce2519 832
a84f14ab 833<tr><td valign="top">
834@anchor Walter07 **[Walter07]**
835</td><td>
27ce2519 836Bruce Walter, Stephen Marschner, Hongsong Li, Kenneth Torrance,
837"Microfacet Models for Refraction through Rough Surfaces", *Proceedings of Eurographics Symposium on Rendering*, 2007.
a84f14ab 838</td></tr>
27ce2519 839
a84f14ab 840<tr><td valign="top">
841@anchor Cao15 **[Cao15]**
842</td><td>
27ce2519 843Jiayin Cao, "Sampling microfacet BRDF", November 1, 2015,
844*Blog post*: [https://agraphicsguy.wordpress.com/2015/11/01/sampling-microfacet-brdf/](https://agraphicsguy.wordpress.com/2015/11/01/sampling-microfacet-brdf/).
a84f14ab 845</td></tr>
27ce2519 846
a84f14ab 847<tr><td valign="top">
848@anchor Schutte18 **[Schutte18]**
849</td><td>
27ce2519 850Joe Schutte, "Sampling techniques for GGX with Smith Masking-Shadowing: Part 1", March 7, 2018,
851*Blog post*: [https://schuttejoe.github.io/post/ggximportancesamplingpart1/](https://schuttejoe.github.io/post/ggximportancesamplingpart1/).
a84f14ab 852</td></tr>
27ce2519 853
a84f14ab 854<tr><td valign="top">
855@anchor Colbert07 **[Colbert07]**
856</td><td>
27ce2519 857Mark Colbert, Jaroslav Krivanek, "GPU-Based Importance Sampling", *NVIDIA GPU Gems 3*, Chapter 20, 2007.
a84f14ab 858</td></tr>
27ce2519 859
a84f14ab 860<tr><td valign="top">
861@anchor Guy18 **[Guy18]**
862</td><td>
27ce2519 863Romain Guy, Mathias Agopian, "Physically Based Rendering in Filament", *Part of Google's Filament project documentation*:
864[https://google.github.io/filament/](https://google.github.io/filament/Filament.md.html)
a84f14ab 865</td></tr>
27ce2519 866
a84f14ab 867<tr><td valign="top">
868@anchor Aguilar17 **[Aguilar17]**
869</td><td>
27ce2519 870Orlando Aguilar, "Spherical Harmonics", *Blog post*:
871[http://orlandoaguilar.github.io/sh/spherical/harmonics/irradiance/map/2017/02/12/SphericalHarmonics.html](http://orlandoaguilar.github.io/sh/spherical/harmonics/irradiance/map/2017/02/12/SphericalHarmonics.html)
a84f14ab 872</td></tr>
27ce2519 873
a84f14ab 874<tr><td valign="top">
875@anchor Ramamoorthi01 **[Ramamoorthi01]**
876</td><td>
27ce2519 877Ravi Ramamoorthi, Pat Hanrahan, "An Efficient Representation for Irradiance Environment Maps",
878*SIGGRAPH '01: Proceedings of the 28th annual conference on Computer graphics and interactive techniques*, August 2001, pp. 497-500
a84f14ab 879</td></tr>
880
881</table>