Introduction
The first game that made with the iPhone game engine described in this series is released before Christmas on App Store, it is called Monster Bowling.

You may surprise that the final product is different from the screen shots shown before which is an adventure game that controlling a ship to explore the world. This is because someone in that project lost his passion and the production of that game has been stopped. And I do not want to waste my effort on making an engine without producing any game, so I decided to work with another artist which resulted in Monster Bowling. The engine and tools took me 9 months to write and the game code took me another 6 months of spare time to complete, which is quite a long time...

What we did wrong
Here I will summarize what mistake we have taken in these 2 projects.

First, the ship game project failed because the scale of the game is too big. Although some of the content is cut during production, there are still too many things need to take care and the team member gradually lose their passion about the game. So, my advice of making a game in spare time, especially when you work with other people the first time, is don't try to make a big game. If anyone in the team lose their passion, the game will never get shipped no matter how hard you work.

Second, the other reason why the ship game failed because there is no a clear direction in the team, how the game project should progress and a clear schedule and check list what should be done. Without a direction, the team member will easily get lost and lost their passion.

Third, this point is about the released game Monster Bowling. When the game was just released, the sales was not good. It is because within this team, there is only 1 artist and 1 programmer, both were putting most of their effort on developing the game and didn't pay much attention to marketing the game and let the public know about the game. So when the game was just released, there didn't have any game review to boost the sales up.

Fourth, this point is about technical issue within the engine. Within the engine, both the game runtime and the editors use the same file format which is a binary format. But doing in this way resulted in a problem which cannot merge the files, for example, when the game objects within a game level are modified by 2 people, we are not able to merge the changes. Luckily, as our team is small, we just do not edit the same file at the same time which does not cause a great trouble in the project.

What we did right
After talking about the mistakes, I will talk about some decisions that I think is made correctly.

First, during the production of Monster Bowling, there are constant communication within the team. We talked about the progress of the game every 1 or 2 weeks. This is very important for those who make games in their free time because through talking about the game, we can keep the passion and let other team member know what are going on in the project.

Second, the decision of making editors and Mac version of the game is important. This enable artist to build and test the game level without any programmer interaction. Although the editors are not user friendly, the artist would suggested to add the minimal set of features (e.g. adding zoom extend to objects in the 3D viewport) that he can live with so that he can work in a much faster way. Explicitly saying what features are lacking will help the programmer to improve the tools.

Level editor used to build the game

Third, embedding a scripting language to the game engine is useful. Writing game play code in Lua is much faster than using C++. Also modifying the script file is much faster than compiling C++ code which can have a faster iteration time. Also scripting language can draw a clear separation between engine code and game play code. Thanks to Lua, when migrating from the ship project to Monster Bowling, there does not have much to change in the engine code.

Test level of Monster Bowling which built using Lua

Fourth, as the engine is written from scratch starting from memory management, rendering code, animation. I still use some open source library: Bullet physics library, Lua and SOIL (for decompressing PNG texture). I use them because my goal is to learn how an engine should be architected to interact with different modules. Also I am not strong at physics and scripting, writing my own physics library and scripting language will take a long time and the code would not be as fast as Bullet and Lua. So I choose to use some open source library which have source code so that it is easier to debug.

Conclusion
This iPhone game engine project took a total of 1 year and 3 months of spare time to complete and produced 1 game. As there is only 1 programmer, Some open source library are used and I decided to sacrifice my favorite graphics programming(the engine doesn't perform lighting and using texture color only) to keep the work load low. I think that it is ridiculous to write an engine without making any game, so I worked on Monster Bowling using this engine, but there are some effort is wasted (i.e. world level streaming) when migrating from the ship project. So if someone just want to focus on making games only, I would suggest they license a commercial engine such as Unity or UDK rather than writing their own engine which will cost a large amount of time. But for me, writing my own engine is to learn and most important is for fun.

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:

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)

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 f0is 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.

IntroductionSpherical 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 cwhich 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.

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.

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.

Dual Quaternion as Rigid Transform
Unit Dual Quaternion can represent a rigid transformation (i.e. rotation & translation) like matrix. Like ordinary quaternion, dual quaternion can represent a rotation using ordinary quaternion with the dual part equals to zero:

To representation a translation (tx, ty, tz), the following dual quaternion can be used:

Then, we can combine the above 2 transform into 1 dual quaternion (which is also unit dual quaternion) to represent a rotation followed by a translation

With the above arithmetic operations, we can transform a point p(px, py, pz), using unit dual quaternion like ordinary unit quaternion:

Dual Quaternion - Matrix Conversion
Sometimes it is convenient to have methods to convert between Unit Dual Quaternion and Matrix. Assume we have ordinary quaternion - matrix conversion functions.

Converting Dual Quaternion to Matrix:
Given any unit dual quaternion which is composed of a rotation followed by a translation, we have:

Therefore, we can find the rotation matrix by considering the real part of the dual quaternion while the translation (i.e. tx, ty, tz) can be solved by using system of linear equations by equating coefficients.

Converting Matrix to Dual Quaternion:
Given a rigid transform matrix, we can decompose the rotation matrix into real quaternion as usual and the translation dual quaternion can also be obtained by:
then multiply the translation dual quaternion with rotation dual quaternion will give the answer.

Blending Rigid Transform using Dual Quaternion
Using dual quaternion to represent transform is better when blending multiple transformations, which can be applied to skinning a mesh. From that paper, it suggests using an fast approximation for blending called Dual quaternion Linear Blending (DLB):

note that the norm in the denominator is a dual number, which can be treated as 1 divided by the norm and gives another dual number (using dual number division) so that it can be multiply by a dual quaternion.

where the real number a is called real part and the real number b is called dual part.

Arithmetic operations
Dual number can perform the arithmetic operations as below:

Addition:

Multiplication:

Division:

Finding derivative using Dual Number
The interesting part of dual number is when it is applied to Taylor Series. When substituting a dual number into a differentiable function using the Taylor Series:

This gives a very nice property that we can find the first derivative, f'(a), by consider the dual part of f(a+bε), which can be evaluated using dual number arithmetic.
For example, given a function

we want to find the first derivative of f(x) at x = 2, i.e. f'(2). We can find it by using dual number arithmetic where f'(2) will equals to the dual part of f(2+ε) according to Taylor Series.

Therefore, f'(2)= 8/9, you can verify this by finding f'(x) and substitute 2 into it, which will give the same answer.

Conclusion
By using dual number, we can find the derivative of a function using dual arithmetic. Hence, we can also find the tangent to an arbitrary point, p, on a given parametric curve which is equals to the normalized dual part of the point p. For those who are interested in finding out more about dual number, I recommend to read the presentation "Dual Numbers: Simple Math, Easy C++ Coding, and Lots of Tricks" by Gino van den Bergen in GDC Europe 2009.