9db2bb3aec99f8745d570e5268b4da54d5830cd9
[occt.git] / dox / dev_guides / visualization / pbr_math.md
1 PBR math (rasterization) {#occt_dev_guides__pbr_math}
2 ========================
3 @tableofcontents
4
5 # Preface
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.
8 Before programmable pipeline has been introduced, graphics cards implemented Gouraud shading as part of fixed-function Transformation & Lighting (T&L) hardware blocks.
9 Nowadays, 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
11 PBR models try to fit surface shading formulas into constrains of physical laws of light propagation / absorption / reflection - hence, called "physically-based".
12 There are two main categories of PBR illumination:
13
14  1. Non-real-time renderer (cinematic).
15  2. Real-time renderer.
16
17 The main objective of cinematic renderer is uncompromised quality, so that it relies on ray-tracing (path-tracing) rendering pipeline.
18 Although 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.
21 The main objective of real-time PBR renderer is to be fast enough even on low-end graphics hardware.
22 So 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
24 OCCT 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.
25 This article describes the math underneath PBR shading in OCCT 3D Viewer and its GLSL programs.
26 However, this article does not clarifies related high-level APIs nor PBR material creation pipelines, as this is another topic.
27
28 # Notation
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
49 # Illumination model
50
51 The main goal of illumination model is to calculate outgoing light radiance \f$L_o\f$ along the certain direction.
52 The 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.
53 Considering 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
57 Where \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$.
58 Opaqueness of the surface mentioned earlier is important because in that case hemisphere is enough.
59 More 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.
61 This is mainly due to geometric laws. The rest part of integral is the key of the whole illumination model.
62 BRDF defines it's complexity and optical properties of material.
63 It has to model all light and material interactions and also has to satisfy some following criteria in order to be physical correct:
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
68 It 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
72 Where \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).
73 So 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
77 PBR theory is based on **Cook-Torrance specular BRDF**. 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.
78 If 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.
79 In 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.
80 In 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.
81 Going 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
85 Three parts presented in nominator have its own meaning but can have different implementation with various levels of complexity and physical accuracy.
86 In that paper only one certain implementation is used. The \f$D\f$ component is responsible for **micro faces normals distribution**.
87 It is the main instrument that controls reflection's shape and strength according to **roughness** \f$r\f$ parameter.
88 The implementation with good visual results is **Trowbridge-Reitz GGX** approach used in Disney's RenderMan and Unreal Engine:
89
90 \f[D=\frac{\alpha^2}{\pi(\cos^2\theta_h(\alpha^2-1) + 1)^2}\f]
91
92 Where \f$\alpha = r^2\f$. This square in needed only for smoother roughness parameter control.
93 Without it the visual appearance of surface becomes rough too quickly during the parameter's increasing.
94
95 The second \f$G\f$ component is called **geometric shadowing** or attenuation factor.
96 The point is that micro surfaces form kind of terrain and can cast shadows over each other especially on extreme viewing angles.
97 **Schlick's model** has been chosen as implementation:
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
101 Where \f$k=\frac{\alpha}{2}\f$, which means \f$k=\frac{r^2}{2}\f$ in terms of this paper.
102 But \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.
103 One of this modification will be described later in following chapters.
104
105 The last component \f$F\f$ shows **how much light is reflected from surface** and is called **Fresnel's factor**.
106 The rest amount of radiance might be absorbed or refracted by material.
107 The most accurate expression of it is pretty complicate for calculation so that there is a variety of approximations.
108 The good one with less computation efforts is **Schlick's implementation**:
109
110 \f[F=F_0+(1-F_0)(1-\cos\theta_{vh})^5\f]
111
112 Here \f$F_0\f$ is material's response coefficient at normal incidence (zero angle).
113 Fresnel'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.
114 In 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.
116 The reference value of \f$IOR\f$ for plastic is **1.5**, and this value can be considered as default for all unknown dielectrics.
117 In practice this parameter controls reflectance ability of material.
118 Also 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.
119 This formula might be further propagated onto metals by using \f$F_0\f$ measured specifically for certain metal.
120 It can be considered as some kind of a 'color' of metal and can be stored as albedo parameter \f$c\f$.
121 And 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
125 For pure dielectrics with \f$m=0\f$ exactly Schlick's approximation will be used.
126 For pure metals with \f$m=1\f$ it will be a little inaccurate but the same formula with measured \f$F_0\f$ values.
127 Everything 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.
128 Intermediate values may represent mixed areas for smooth transition between materials - like partially rusted metal (rust is mostly dielectric).
129 Also it might be useful when parameters are read from textures with filtering and smoothing.
130
131 BRDF described above has one important trait making computations easier called **isotropy**.
132 Isotropy in this case means independence from rotation about normal resulting from supposition of uniform micro faces distribution at any direction along a surface.
133 It 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.
134 Of course, isotropic materials form only subset of all real world's materials, but this subset covers majority of cases.
135 There are special models considering special anisotropic traits of surfaces like a grinding of metal or other with dependency on rotation about normal;
136 these models require special calculation tricks and additional parameters and are out of scope of this paper.
137
138 The only thing left to do is to define \f$f_d(v,l)\f$.
139 This part is responsible for processes happening in depth of material.
140 First of all the amount of input light radiance participating in these processes is needed to be calculated.
141 And 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
145 This part of ingoing light is assumed to be refracted in depth of surface and variety of events may happen there.
146 A sequence of absorptions, reflections and reemissions more or less leads to light's subsurface scattering.
147 Some part of this scattered light can go back outside but in modified form and in pretty unpredictable directions and positions.
148 For opaque materials this part is noticeable and forms it's own color.
149 If 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.
150 This assumption covers a big amount of real world opaque materials.
151 Other materials like skin, milk etc. with noticeable effect of subsurface scattering usually presented in form of partial translucency and some kind of self emission
152 have more widely distributed output points and require more accurate and complicate ways of modeling with maybe some theory and techniques from volumetric rendering.
153 The 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
157 Where \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.
158 There is one detail about light interaction bringing some physicality to the model, and that is an absence of this diffuse component in metals.
159 Metals reflect main part of light and the rest of it is absorbed being transformed into other form (mostly heat).
160 That 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
162 So 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
168 In this chapter one possible implementation of illumination model reflecting main PBR principles has been defined.
169 The next step is using of it in practice.
170
171 # Practical application
172
173 It's time to apply deduced illumination model in practice.
174 And the first step of it is separation of **direction based light sources** from illumination integral.
175 Directional nature of such light sources means possibility to calculate it's influence to point of surface using only one direction and its intensity.
176 Usually 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).
177 This is just a kind of abstraction, while real world light emitters have noticeably sizes.
178 But 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.
179 In 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.
180 Having 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
184 Where \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.
186 The BRDF can be used directly without any calculation problems.
187 The only exception might be \f$k\f$ in \f$G\f$ factor - it is recommended to be equal \f$\frac{(r+1)^2}{8}\f$ in order to get more pleasant results (that is modification mentioned in previous chapter).
188 And actually it is enough to finally see something.
189 There will be correct visualization with assumption of complete dark environment and absence of other points influence.
190 It 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
194 It includes influence of light reflected or scattered from other points and environment's contribution.
195 It's impossible to achieve photorealistic results without this component, but is is also very difficult to compute.
196 While 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.
197 But right now lets summarize the practical application of illumination model.
198 At this moment the output radiance is represented as:
199
200 \f[L_o = L_{direct} + L_{indirect}\f]
201
202 Where \f$L_{direct}\f$ is direction based light sources contribution which can be directly computed just applying bare formulas.
203 It 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.
205 But 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
207 # Image based lighting
208
209 The next goal after \f$L_{direct}\f$ calculation is to find \f$L_{indirect}\f$.
210 And it would be easier if \f$L_i^{indirect}(l)\f$ was known for every \f$l\f$.
211 That is the main assumption of **image based lightning** (**IBL**).
212 In practice, it can be achieved using environment image map, which is a special image representing illumination from all possible directions.
213 This 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.
214 Environment image might be packed in different ways - **cube maps** and equirectangular maps are the most commonly used.
215 Anyway, 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.
216 Lets 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
220 Substituting 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
226 This splitting seems not to lead to simplicity of calculation but these two parts will be computed in slightly different ways in future.
227 Lets 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
233 Next transformations of these expressions require understanding of numerical way to find hemisphere integral and also its performance optimization techniques.
234 And that the topic of the next chapter.
235
236 # Monte-Carlo numeric integration
237
238 **Monte-Carlo** is one of numeric methods to **find integral**.
239 It is based on idea of mathematical expectation calculation.
240 In 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
244 Also 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
248 It can be used in general definite integrals calculations.
249 Just 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
253 Where \f$f(x)\f$ is considered to be equal to zero outside of \f$[a, b]\f$ range.
254 This 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
258 The main questions are choosing \f$p(l)\f$ and generating samples \f$l_i\f$.
259 The one of the simple ways is uniform distribution along sphere or hemisphere.
260 Lets realize that on sphere for example.
261 There 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
265 Where \f$S\f$ is the unit sphere.
266 In 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
270 So that \f$p(l) = \frac{1}{4\pi}\f$.
271 Usually 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.
272 But 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$.
274 The 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
278 Where \f$r = 1\f$.
279 In 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
283 Where:
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
287 So that joint probability density in new variables looks like:
288
289 \f[p(\theta, \phi) = \frac{\sin\theta}{4\pi}\f]
290
291 This variable transfer rule of **Probability Density Function** (**PDF**) will be useful in following chapters, when integral calculation optimization techniques will be being told.
292 Having \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
298 The final step is sequence generation itself.
299 In 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$.
300 And that can be done via any known true- and pseudo- random generators.
301 Even simple \f$\frac{1}{i}\f$ sequence is appropriate for beginning but it can be not so efficient in terms of computations convergence.
302 There are specially designed series for the last reason and it will be tackled in chapter about optimizations.
303 The \f$\phi\f$ variable is noticed to be uniformly distributed so that it can be directly generated without any additional manipulations.
304 Just range \f$[0, 1]\f$ is needed to be mapped to range \f$[0, 2\pi]\f$.
305 For any other variables including \f$\theta\f$ the inverse transform sampling approach can be applied.
306 First of all **cumulative distribution function** (**CDF**) is needed to be found.
307 It is probability of random value to be less than argument of this functions by definition.
308 For 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
312 Lets 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
316 The CDF maps \f$\theta\f$ values from range of \f$[0, \pi]\f$ to probability in range of \f$[0, 1]\f$.
317 The 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
321 If 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.
322 In 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
326 That is the key of this random values generation technique.
327 All 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
339 Mote-Carlo integration cannot guarantee exact estimation of convergence speed with using random generated samples.
340 There is only probability estimation of it.
341 But this algorithm is pretty universal and relatively simple to be applied to almost any function using at least uniform distributed points.
342 Moreover 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).
343 That is why this method is widely used in computer graphics and demonstrates good results.
344 Also 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
346 # Split sum
347
348 Lets go back to the image based lighting and the figure of specular component.
349 As 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
353 The 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.
358 Optimization strategies use different samples distributions for different view direction orientations and roughness values.
359 Anyway even with all optimization techniques this algorithm continues to require too much calculations.
360 Good visual results require noticeable number of samples and using this approach for every point in real time rendering becomes unrealistic.
361 The way to avoid these enormous calculations is doing them beforehand somehow.
362 The first trick on the way to this is split the sum separating environment light component:
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
366 Where 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
370 This integral is exact \f$L_{indirect}^s\f$ in condition when \f$L_i^{indirect}(l) = 1\f$ what just means white uniform environment.
371 The sum before it is kind of averaged environment illumination.
372 The main accomplishment after all this manipulations is possibility to calculate light and BRDF components separately.
373 The 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.
374 The 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.
375 Variation of normal is not enough in that case and these variables are needed to be considered too.
376 The 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.
377 And 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.
378 That 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.
379 For example completely mirroring materials with \f$r = 0\f$ will not be looked as expected if just uniform distribution is used
380 because 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
382 The 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
386 This integral is not actually a scalar.
387 That is RGB value due to only \f$F\f$ factor and even more only to \f$F_0\f$.
388 In order to simplify future computations \f$F_0\f$ is needed to be moved out of integral.
389 Substitution 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
395 This form may not look easier, but it has several advantages.
396 The 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).
397 The 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.
398 And 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$.
399 Considering 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
401 Current result for \f$L_{indirect}^s\f$ is computing it using 2D image for BRDF part and omnidirectional image for environment illumination.
402 There 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.
403 It's time to talk about that.
404
405 # Importance sampling
406
407 Current 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
411 Where \f$g(v, l)\f$ is just arbitrary function representing Fresnel's factor itself or its components.
412 In 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.
413 So 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.
414 That 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.
416 And it will be like that if PDF is already part of original function because in that case probability density directly affects to contribution forming.
417 Separating this distribution component is one possible and effective way to realize importance sampling strategy.
418 In integral presented above PDF part already exists and that is \f$D\f$ component.
419 But 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.
420 Anyway 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
425 Frankly speaking \f$D(h)\f$ is called normal distribution but cannot be directly used as hemisphere distribution.
426 Originally it is statistical factor used to define total area of micro faces \f$\mathrm{d}A_h\f$
427 whose normals lie withing 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$:
428
429 \f[dA_h = D(h)\,\mathrm{d}h\, \mathrm{d}A\f]
430
431 First of all this factor must be positive:
432
433 \f[D(h) \geq 0\f]
434
435 But 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
439 This trait does not allow to use \f$D\f$ as hemisphere distribution.
440 But 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
444 Which means that total area of micro faces projected to any direction must be the same as projected area of origin macro surface.
445 It 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
449 So that \f$\cos\theta_h\f$ coefficient normalizes normal distribution in terms of hemisphere and allows to use it as distribution.
450 Finally 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.
455 Lets 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.
462 The fact of independency from \f$\phi\f$ is enough to figure out that this coordinate is uniformly distributed.
463 Anyway the \f$F(\theta_h)\f$ is next step:
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
469 In 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
473 So that there is no more obstacles to generate \f$h\f$.
474 But the main goal was \f$l\f$ direction.
475 In 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
479 That is practical side of light direction generation.
480 But the theoretical one is needed to be resolved to calculate sum.
481 It is time to find \f$p(l)\f$ using known \f$p(h)\f$.
482 First of all the fact that \f$l\f$ is just transformed \f$h\f$ is needed to be understood.
483 In that way the light direction's PDF has following expression:
484
485 \f[p(l) = p(h)|J_T|\f]
486
487 Where \f$|J_T|\f$ is Jacobian of reflection transformation.
488 Lets find it.
489 Right now \f$n\f$ is axis from which \f$\theta\f$ spherical coordinate is encountered.
490 The first step is setting \f$v\f$ as starting point of \f$\theta\f$ instead of \f$n\f$.
491 This is linear transform so that \f$|J_T| = 1\f$.
492 Next step is transfer to spherical coordinate with \f$|J_T| = \sin\theta_{vh}\f$.
493 Due to previous step \f$\theta_{vh}\f$ is used instead of \f$\theta_h\f$.
494 In 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$.
495 In series of transform the Jacobians are multiplied so that currently \f$|J_T| = \frac{1}{2}\sin\theta_{vh}\f$.
496 And the final step is inverse transform to Cartesian coordinate system with \f$|J_T| = (\sin\theta_{vl})^{-1} = (\sin2\theta_{vh})^{-1}\f$.
497 Combining this all together the following expression is obtained for reflection transform Jacobian:
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
501 And 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
505 Lets 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
509 Here \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.
510 Of 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
513 Summarizing importance sampling strategy described above the convergence of Monte-Carlo integration can be improved using special PDF correlated with integrated function.
514 In case of BRDF with normal distribution functions \f$D\f$ the PDF producing procedure is defined.
515 Practically half vector \f$h\f$ is generated first and \f$l\f$ is obtained from it by view vector \f$v\f$ reflecting.
516 Due to this transformation final form of probability density used in sum is quite different but also has defined algorithm of calculation.
517 For 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.
518 The same samples generation procedure must be used in specular map baking described in next chapter.
519
520 # Specular map
521
522 The 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.
523 That 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
527 As was mentioned this sum must be calculated for every normal direction using the same samples generation principles as in numeric integration computation.
528 This principles require two scalar parameters \f$\cos\theta_v\f$ and \f$r\f$ but now \f$\phi\f$ really matters.
529 So that in fact the specular map has to be saved in 3D table consisting omnidirectional textures.
530 That is a big expense of computational and memory resources.
531 A couple of tricks helps to reduce dimensions.
532 First of all the \f$\cos\theta_v\f$ and \f$\phi\f$ can be just excluded.
533 In that way \f$v\f$ is considered to be equal to \f$n\f$.
534 Of course this approach produces an error and affects the final result.
535 It can be fixed more or less by \f$\cos\theta_{l_i}\f$ weighting:
536
537 \f[\frac{1}{N} \sum_{i=1}^N L_i^{indirect}(l_i) \cos\theta_{l_i}\f]
538
539 It 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.
540 The problem is that for \f$v \neq n\f$ considering this sum to be defined related to \f$n\f$ became incorrect.
541 For example, for complete mirroring materials with \f$r = 0\f$ this sum must boil down to \f$L_i^{indirect}(v_r)\f$
542 but 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$.
543 That it just mirroring reflection principle.
544 Assumption of \f$n = v\f$ also means that \f$n = v = v_r\f$.
545 In that way radiance map better to be considered as averaging of illumination coming from \f$v_r\f$.
546 So that it has become to be defined related to reflection direction which has to be calculated before map's fetching.
547
548 Anyway, there are just two dimensions in radiance look-up table remain.
549 The rest one with \f$r\f$ parameter cannot be reduced.
550 There 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.
551 This is another source of visual artifacts but it also works good in practice and that is pretty common approach.
552 But it still requires noticeably amount of samples and that is for every pixel related to each \f$r\f$ value.
553 It 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.
554 And there are several possible ways to improve radiance map baking performance.
555
556 The first one is using textures with smaller resolutions for larger roughnesses.
557 The 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.
558 So less pixels in combination with linear interpolation is enough to store less detailed convolutions.
559 Moreover, this approach naturally works with textures levels of details in graphics API
560 so 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.
561 As practice shows 6 levels are enough.
562
563 After reducing pixels count it is turn for samples number.
564 And again correlation with roughness can be noticed.
565 For example map for completely mirroring materials with \f$r = 0\f$ the same sample \f$l_i = v_r\f$ will be produced.
566 So that only one sample is enough in this case.
567 In opposite way directions for \f$r = 1\f$ will be scattered over almost whole hemisphere what requires as much samples as available.
568 The '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.
569 But 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.
571 In 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$.
572 The inverse task of range searching using given \f$u\f$ can be solved with help of \f$F^{-1}(u) = \theta'_h\f$.
573 If \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$)
574 then \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.
575 After 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
579 It 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).
581 This pretty tricky technique allows reducing calculations for every pixels without sacrificing the quality.
582
583 In addition to optimizations mentioned before another one can be applied in order to help to reduce samples numbers as previous one.
584 Using less samples produces image noise due to discrete nature of Monte-Carlo approximation.
585 But it can be slightly smoothed using some prefiltration.
586 The 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.
587 So that this direction represents the averaged illumination from relatively big area on hemisphere but approximate it by just a constant.
588 It wold be better to get from such direction already averaged over bigger area environment.
589 It 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.
590 Also mip levels generation is build in most common graphic API so there are no problems with it.
591 But first of all the area covered by one sample is needed to be found.
592 And that can be done as:
593
594 \f[\Omega_s = \frac{1}{N\,p(l)} = \frac{4\cos\theta_{vh}}{ND\cos\theta_h}\f]
595
596 Circumstance 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
600 Of course all zero divisions must be avoided by clamping, for example.
601 After that the area covered by one pixel of environment map is calculated.
602 In case of a cube map it looks like:
603
604 \f[\Omega_p = \frac{4\pi}{6k^2}\f]
605
606 Where \f$k\f$ is size of one cube map side in pixels (sides are assumed to be quads).
607 Finally 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
611 The mathematics connected with mip levels sizes lie behind it but this is out of scope of this paper.
612 In combination with previous optimization technique this approach allows \f$N(1)\f$ to be smaller keeping visual results good.
613
614 That 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
616 In that way \f$L_{indirect}^s\f$ can be completely computed without any integral approximations.
617 Only 2D look-up table of BRDF part and mip mapped omnidirectional texture of irradiance map are needed.
618 The first one can be got even without any environment.
619 It 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.
620 In order to complete whole image based lighting the \f$L_{indirect}^d\f$ component is left to be discussed.
621
622 # Spherical harmonics
623
624 Lets 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
628 Of course, Monte-Carlo algorithm can be applied directly and hemisphere integral can be precalculated for every normal direction
629 but 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).
630 In order to resolve it modified version of Schlick's approximation has been created:
631
632 \f[F \approx F_{ss}=F_0+(\max(1-r, F_0))(1-\cos\theta_v)^5\f]
633
634 It differs from origin one and loses accuracy a little bit but now there is no light direction inside
635 so 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
639 The resulted expression without \f$(1-m)\f$ and \f$(1-F_{ss})\f$ parts is pretty known entity called **irradiance**.
640 It 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).
641 But even with that samples will be scattered almost over whole hemisphere.
642 As was discussed in previous chapter this case requires significant amount of samples in order to average illumination with appropriate quality.
643 Poor accuracy resulted from lack of summand can be noticed especially on high frequency environments having a lot of contrasting details.
644 It worth to be mentioned that irradiance is used not only in BRDF.
645 Omnidirectional diffuse illumination captured for certain point or even for several points uniformly or hierarchically distributed is base of some baking global illumination techniques.
646 There it is called a **light probe**. So that other way to compute and store irradiance maps was found resolving many mentioned problems.
647 The Fourier's decomposition analogue for spherical function allows to achieve this.
648 That would be easy to explain concept directly on example.
649 So lets start from \f$L_i^{indirect}(l)\f$.
650 It is spherical function because directions are just points on sphere.
651 The 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
655 Where \f$y_i^j(l)\f$ are spherical functions forming orthonormalized basis called spherical harmonics and \f$L_i^j\f$ is decompositions coefficients.
656 Orthonormality means that dot product of two basis elements is equal to 1 if this is the same functions and is equal to zero otherwise.
657 Dot 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
661 Function basis with such traits is known and is described by following formulas:
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
673 Here \f$K_i^j\f$ are normalization factors and \f$P_i^j\f$ are **Legendre's polynomials**.
674 Decomposition 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
678 Fact that all calculation happen over a sphere but not over hemisphere right now is important not to be missed.
679 That 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
683 First 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
687 Where \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
691 Resulted 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.
692 This 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.
693 That is pretty tricky and explanation about meaning of convolution on sphere is out of scope of this paper.
694 Fact that this is convolution analogue related to \f$n\f$ is enough for now.
695 The 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
700 Where \f$c_i^0\f$ are spherical harmonics factors corresponding to \f$\overline{\cos}\theta\f$.
701 Therefore 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.
706 Due to independence from \f$\phi\f$ all \f$c_i^{j \neq 0} = 0\f$.
707 The 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$.
714 There 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
720 Starting 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.
721 The same principle is true for convolution too because its coefficients are multiplied by \f$c_i\f$.
722 So there is no need to use more than even three bands (\f$i = 0, 1, 2\f$) of basis functions.
723 Lets write down them all in Cartesian coordinate system:
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
737 All \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.
738 Other 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
742 Finally 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
746 Where \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.
747 Of 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.
748 Moreover, texture is not needed at all in that case and only 9 colors representing \f$L_i^j\f$ have to be saved.
749 The 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.
751 So there is sense to store only squares of it.
752 All tables with constants presented below:
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
768 Summarizing all mathematics above spherical harmonics decomposition boils down irradiance map to only 9 values which is needed to be precalculated as integrals.
769 As practice shows this is very good approximation of diffuse indirect illumination component.
770
771 # Transparent materials
772
773 TODO
774
775 # Low discrepancy sequence
776
777 TODO