Extracting dominant light from Spherical Harmonics

Introduction
Spherical Harmonics(SH) functions can represent low frequency data such as diffuse lighting, where those high frequency details are lost after projected to SH. Luckily we can extract a dominant directional light from SH coefficients to fake specular lighting. We can also extract more than 1 directional light from SH coefficients, but this post will only focus on extracting 1 dominant light, those interested can read Stupid Spherical Harmonics (SH) Tricks for the details. A webGL demo is provided at the last section which will only extract 1 directional light.

Extracting dominant light direction
We can get a single dominant light direction from the SH projected environment lighting, Le. Consider we approximate the environment light up to band 1 (i.e. l=1):

Finding the dominant light direction is equivalent to choose an incoming direction, ω, so that Le(ω)is maximized. In other words, cosθ should equals to 1:

So we can extract the dominant light direction for a single color channel. Finally the dominant light direction can be calculated by scaling each dominant direction for RGB channels using the ration that convert color to gray scale:

Extracting dominant light intensity
After extracting the light direction, the remaining problem is to calculate the light intensity. That's mean we want to calculate an intensity s, so that the error between the extracted light and the light environment is at minimum (Le is the original environment light while Ld is the directional light):

To minimize the error, differentiate the equation and solve it equals to zero:

If both lighting functions are projected into SH, the intensity can be simplified to:

The next step is to project the directional light(with unit intensity) into SH basis (ci is the SH coefficient of the projected directional light):

Therefore the SH coefficients of projected directional light can be calculated by substituting the light direction into the corresponding SH basis function.

As the SH projected directional light is in unit intensity, we want to scale it with a factor so that the extracted light intensity s is the light color that can be ready for use in direct lighting equation which is defined as (detail explanation can be found in [4]):

For artist convenience, clight does not correspond to a direct radiometric measure of the light’s intensity; it is specified as the color a white Lambertian surface would have when illuminated by the light from a direction parallel to the surface normal (lc = n).

So we need to calculate a scaling factor, c, that scale the SH projected directional light such that:

We can project both L(ω) and (n . ω) into SH to calculate the integral. To project the transfer function (n . ω) into SH, we can first align the n to +Z-axis, which is zonal harmonics, then we can rotate the ZH coefficient into any direction using the equation:

The ZH coefficients of (n . ω) are: (note that the result is different from Stupid Spherical Harmonics (SH) Tricks in the Normalization section as we have taken the π term outside the integral)

Then rotate the ZH coefficients such that the normal direction is equals to the light direction, ld (because we need ld = n as stated above), we have:

Finally we can go back to compute the scaling factor, c,  for the SH projected directional light (we calculate up to band=2):

Therefore the steps to extract the dominant light intensity are first to project the directional light into SH with a scaling factor c, and then light color, s,  can be calculated by:

WebGL Demo
A webGL demo (need a webGL enabled browser such as Chrome) is provided to illustrate how to extract a single directional light to fake the specular lighting from the SH coefficient. The specular lighting is calculated using the basic Blinn-Phong specular team for simplicity reason, other specular lighting equation can be used such as those physically plausible. (The source code can be downloaded from here.)

Your browser does not support the canvas tag/WebGL. This is a static example of what would be seen.

Render Diffuse
Render Specular
Rotate Model
Glossiness

Render Diffuse

Render Specular

Rotate Model

Glossiness

Conclusion
Extracting the dominant directional light from SH projected light is easy to compute with the following steps: First, calculate the dominant light direction. Second, project the dominant light into SH with a normalization factor. Third, calculate the light color. The extracted light can be used for specular lighting to give an impression of high frequency lighting.

References
[1] Stupid Spherical Harmonics (SH) Tricks: http://www.ppsloan.org/publications/StupidSH36.pdf

[2] Light Factorization for Mixed-Frequency Shadows in Augmented Reality: http://zurich.disneyresearch.com/~wjarosz/publications/nowrouzezahrai11light.pdf
[3] Physically Based Rendering, Second Edition: From Theory To Implementation Ch.17.2.2
[4] Physically-Based Shading Models in Film and Game Production: http://renderwonk.com/publications/s2010-shading-course/hoffman/s2010_physically_based_shading_hoffman_a_notes.pdf

[5] PI or not to PI in game lighting equation: http://seblagarde.wordpress.com/2012/01/08/pi-or-not-to-pi-in-game-lighting-equation/
[6] March of the Froblins: Simulation and Rendering Massive Crowds of Intelligent and Detailed Creatures on GPU: http://developer.amd.com/documentation/presentations/legacy/Chapter03-SBOT-March_of_The_Froblins.pdf
[7] Pick dominant light from sh coeffs: http://sourceforge.net/mailarchive/message.php?msg_id=28778827

Microfacet BRDF

Introduction
In recent years, there are more and more papers talking about applying physically based BRDF in games. So I decided to spend some time to investigate it. For a BRDF to be physically plausible, it should satisfy 2 conditions:

1. Reciprocity: The incident light direction(l) and reflected light direction(r) for a BRDF(f) is the same after the incident and reflected direction is swapped. i.e. f(l, r)= f(r, l)
2. Energy Conservation: The total energy of reflected light is less than or equal to the energy of the incident light. i.e.

                             

A physically based specular BRDF is based on micro-facet theory, which describe a surface is composed of many micro-facets and each micro-facet will only reflect light in a single direction according to their normal(m):

So, in the above diagram, for light coming from direction l to be reflected to viewing direction v, the micro-facet normal m must be equals to the half vector between l and v.
A micro-facet BRDF has the following form:

which consists of 3 terms: Fresnel term(F), Distribution term(D) and Geometry term(G). Their meaning can be found in the background talk presented by Naty Hoffman in siggraph 2010. And these 3 terms can be chosen independently as stated in the talk Physically-based lighting in Call of Duty:Black Ops (although "Microfacet Models for Refraction through Rough Surfaces" states that some G depends on D to maintain energy conservation, but some G are extended to handle arbitrary distribution, so in this blog post, I assume that the G function is independent of D). So I decided to find some distribution functions D and geometry functions G and play with different combinations to see how it affects the rendering result. You can also play around with different combinations using the WebGL demo(need a webGL enabled browser such as Chrome) in the last section of the post.

Fresnel Term
In this test, I use the common Schlick approximation to the Fresnel equation:

and f0 is found by using the following equation:

where the refractive index n can be tuned in the demo.

Distribution Term
Distribution term is used to describe how the microfacet normal distributed around a given direction. In the demo, I used two distribution function: Blinn-Phong and Beckmann distribution function.

For Blinn-Phong distribution, we can derive the distribution function by satisfying the equation:

which means that the projected microfacet area is equal to macro surface area for any projected direction v. So we choose v=n which simplify the equation:

To derive the Blinn Phong distribution function from original Blinn Phong specular term, we just need to multiply a constant K to satisfy the equation:

While the Beckmann distribution has the following form:

To convert between the roughness m in Beckmann distribution and shininess α in Blinn-Phong distribution, the following formula is used:

which gives a very similar result when both refractive index n and roughness m are small. When n>10 and m>0.5 the 2 distribution start to show difference and the difference will get larger when both m and n are getter larger.

Geometry Term
Geometry term is used for describing how much the microfacet is blocked by other microfacet. In the demo, 4 geometry terms have been tested: implicit, Cook-Torrance, Schlick approximation to Smith's shadowing function and Walter approximation to Smith's shadowing function.

The first one is implicit geometry function which has the form:

It is called implicit because when it is used, the microfacet BRDF will only depends on Fresnel equation and distribution function.

The second one used for testing is Cook-Torrance geometry function:

And the other 2 geometry functions used are both trying to approximate the Smith's shadowing function which decompose the geometry function into another 2 geometry function as below:

With Schlick's approximation, the following G1 is used:

While Walter's approximate G1 as:

Among 4 geometry terms, the implicit one always show a darker specular color. While the other 3 geometry functions have similar appearance when the roughness m is small. When m is getter larger, the Schlick function will slightly darker than the Cook-Torrance and Walter geometry function. Both Cook-Torrance and Walter function gives a very similar results:

Energy Conservation between Diffuse and Specular BRDF
Energy conservation is important for a physically based BRDF, but most paper only talks about the conservation within the specular BRDF. How about the energy conservation between the diffuse term and specular term? I can only find 2 ways to do this from the paper provided by Tri-Ace. They multiply the diffuse reflection term with a diffuse Fresnel term:

And they later discovered that this term can be approximated with (1- f0), which will show very similar results.  However, using this term will violate the reciprocity of the BRDF. If diffuse energy conservation is enabled, when the refractive index change, the ratio between the diffuse and specular reflection also change accordingly.

WebGL Demo
I provide a webGL program so that you can play around with the settings I described above. The model is illuminated by a single white directional light and the red color is the diffuse color. The diffuse BRDF is just a lambert surface which can be turned off in the demo. Dragging inside the viewport can rotate the camera. The source code can be downloaded from here.

Your browser does not support the canvas tag/WebGL. This is a static example of what would be seen.

Distribution Term: Beckmann Blinn Phong	Geometry Term: Implicit Cook Torrance Schlick Walter	Diffuse Energy Conserved: None 1-Fdiff 1-f0
Refractive Index, n (1.0 - 20.0)	(n= / Fresnel 0: )
Roughness, m (0.01- 1.0)	(m=)
	Render Diffuse	Rotate Model

Conclusion
Physically plausible BRDF can give a different material appearance for a surface compare to traditional lighting model. However, in this post, I only use 1 microfacet BRDF for all 3 RGB channels, using different BRDF settings for difference channels is also possible as some material like copper and gold have different f0 term in RGB channels. Also only direct lighting is investigated where secondary lighting BRDF will be left for future blog post.

Reference
[1] Background: Physically-Based Shading (Naty Hoffman): http://renderwonk.com/publications/s2010-shading-course/hoffman/s2010_physically_based_shading_hoffman_a_notes.pdf
[2] Practical Implementation of Physically-Based Shading Models at tri-Ace (Yoshiharu Gotanda): http://renderwonk.com/publications/s2010-shading-course/gotanda/course_note_practical_implementation_at_triace.pdf
[3] Crafting Physically Motivated Shading Models for Game Development (Naty Hoffman): http://renderwonk.com/publications/s2010-shading-course/hoffman/s2010_physically_based_shading_hoffman_b_notes.pdf
[4] Physically-based lighting in Call of Duty: Black Ops: http://advances.realtimerendering.com/s2011/Lazarov-Physically-Based-Lighting-in-Black-Ops%20(Siggraph%202011%20Advances%20in%20Real-Time%20Rendering%20Course).pptx
[5] Microfacet Models for Refraction through Rough Surfaces:http://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.pdf
[6] http://www.rorydriscoll.com/2009/01/25/energy-conservation-in-games/

[7] http://seblagarde.wordpress.com/2011/08/17/hello-world/

Spherical Harmonic Lighting

Introduction
Spherical Harmonics(SH) functions are a set of orthogonal basis functions defined in spherical coordinates using imaginary numbers. In this post, we use the following conversion between spherical and cartesian coordinates:

Since we are dealing with real value functions, we only need to deal with real spherical harmonics functions which in the form of:

The index l of the SH function is called the band index which is an integer >= 0 and index m is an integer with range -l<=m<=l , so there will be (2l + 1) functions in a given band. You may refer to the Appendix A2 of Stupid Spherical Harmonics(SH) Trick to look up the evaluated value of the SH basis function for a pair of (l, m).

The linear combination of SH basis functions with scalar values can be used to approximate a function as below:

With an approximation up to band l = n - 1, which n×n coefficients are needed.
So the remaining problem to approximate a function is to compute the coefficient c which can be solved either analytically or numerically by Monte Carlo Integration.

Monte Carlo Integration

To compute a definite integral numerically, we can consider the Monte Carlo Estimator:

When the number of samples, N, is large enough, the estimator F will equal to the definite integral because considering the expected value of F:

When number of samples,N, is large enough, by the law of large numbers, the estimator F will converge to the definite integral. Therefore, we can calculate the coefficient of the SH basis functions by using Monte Carlo Estimator.

Properties of Spherical Harmonics Function

There are 2 important properties properties of SH functions:

First, it is rotationally invariant. 

Where the rotated function g is still a SH function which its coefficients can be computed by using the coefficients of f. For details of rotating a general SH functions, you can refer to the section 'Rotating Spherical Harmonics' in Spherical Harmonics Lighting: The Gritty Details.

Second, when integrating 2 SH projected functions over the spherical domain, the results will equals to dot product of their SH coefficients (due to the SH basis functions are orthogonal):

This is a nice property that we can calculate the integration over the spherical domain by a dot product of the SH coefficients.

Lighting with SH functions

When performing lighting calculation, we need to solve the rendering equation:

For shading lambert diffuse surface without shadow, we can simplify the rendering equation into:

To solve this integral, we can project the functions L(x, ω) and max(N.ω, 0) into SH functions using Monte Carlo Integration, then by the property 2 described above, the integral equals to dot product of the SH coefficients of the 2 SH projected functions.

Zonal Harmonics
If a SH projected function is rotational symmetric about a fixed axis, it is called Zonal Harmonics(ZH). If this axis is the z-axis, this will make the ZH function only depends on θ, which will result in only one non-zero coefficient in each band with m= 0. Then rotation of the ZH function can be greatly simplified. When the ZH function is rotated to a new axis d, the coefficients of the rotated SH function will equals to:

,which is faster than the general SH rotation. The ZH function is well suit to approximate the function max(N.ω, 0) in the above diffuse surface rendering equation since the SH projected L(x, ω) is usually done in world space while the shading surface can be re-oriented to the same space to perform lighting calculation.

WebGL Demo
Below is a webGL demo (which need a webGL enabled browser such as Chrome) using the cube map on the right as light source and projected to SH function using Monte Carlo Integration.

Both the white and the blue color on the model is reflected from the sun and the blue sky using SH coefficients generated from the cube map and the ZH coefficients projected from max(N.ω, 0) which rotated to world space according the surface normal. The approximation is done up to band l=2.  You can drag in the viewport to rotate the camera.

Your browser does not support the canvas tag/WebGL. This is a static example of what would be seen.

The source code of the webGL can be downloaded here.

Conclusion
SH functions can be used to approximate the rendering equation with only a few coefficients and a simple dot product to evaluate lighting during run time. But it also has its disadvantage while SH can only approximate low frequency function as it needs large number of bands to represent high frequency details.

Reference
[1] Spherical Harmonics Lighting: The Gritty Details: http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf
[2] Stupid Spherical Harmonics(SH) Trick: http://www.ppsloan.org/publications/StupidSH36.pdf
[3] Physically Based Rendering: http://www.amazon.com/gp/product/0123750792/ref=pd_lpo_k2_dp_sr_1?pf_rd_p=486539851&pf_rd_s=lpo-top-stripe-1&pf_rd_t=201&pf_rd_i=012553180X&pf_rd_m=ATVPDKIKX0DER&pf_rd_r=09AG8FQQWKJHC2AEFPD1
[4] Sky box texture downloaded from: http://www.codemonsters.de/home/content.php?show=cubemaps