Solving Global Illumination using the Radiosity Model

 

The radiosity model is defined by a Fredholm Integral equation of the second kind:

This equation represents the stable distribution of light within environment and needs to be solved to solve the global illumination problem. 

These are the steps I followed to solve the equation:

To solve the problem, I used the finite element method to discretize my radiosity model into a sum of easier polynomial functions of the form:

The basis functions 𝑁𝑖(𝑥) I used are constant basis functions which are 1 within element and 0 otherwise.

Galerkin Form of Weighted Residuals.

I needed to define an error metric to approximate my radiosity model, I chose the Galerkin Form of weighted residuals to minimize the error as an average across element surface, so:

Where 𝑊𝑖(𝑥) is a weighting function defined to be the same basis function as 𝑁𝑖(𝑥), so expanding 𝑟(𝑥):

And expanding our approximation function 𝐵̂(𝑥):

Grouping terms:

This is the system of equations to solve. This system of equations can be expressed in matrix form as:

 

The entries of K are given by:

And the entries in E are:

We can simplify the above equation by assuming constant basis functions, and constant reflectivities and emissivities within element.

 

Where 𝛿𝑖𝑗 is the Kronecker Delta function:

 

Then:

Where 𝐸𝑖 is the average area emission value for element 𝑖 . 

Then we can change the integral, to integrate over element areas only: 

Making these substitutions in equation 6:

Dividing through 𝐴𝑖 and moving emission term to the right side gives:

Which is:

The Form-Factor 

The form factor (𝐹𝑖𝑗 equation 16), specifies the fraction of the energy leaving one surface which lands in another. 

The Form Factor between differential areas is defined by:

The Hemi-Cube Algorithm

If the distance between two elements is much larger than the sizes of each element, I can integrate over Area 𝑗 only. For that I used equation 17 and set a hemi-cube around the normal in patch 𝑖 and meshed each of its faces into square pixels (100 x 100 per face).

 

The hemicube contains 5 faces, 1 full face and 4 half-faces. The area of each pixel determines the multiple to multiply the form factor by, so the form factor for each pixel can be calculated as:

 

Each 𝐹𝑝𝑖𝑥𝑒𝑙 value goes into a matrix of form factors. (e.x. if we are on element 𝑖 and we see element 𝑗 through pixel 𝑘 then 𝐹𝑝𝑖𝑥𝑒𝑙 𝑘 value goes into row 𝑖 column 𝑗 of Form Factor Matrix.) 

Making this calculation for all elements in scene, we get filled a full n by n Form Factor Matrix, the resulting matrix can be seen in the next image (NOTE: it is possible that the values cannot be seen because they are relatively low values, but that is good because the sum of each row needs to be less than one to make the system converge): 

NOTE: The matrix images are color coded. Red means negative value, Blue means positive value greater than 1.0 and gray scale means between 0 and 1 with 1 meaning white and 0 meaning black.

The Energy Balance solution

Now having the form-factor matrix we can solve the system of linear equations:

 

Where K is defined to be:

M is the identity matrix; P is a diagonal matrix with reflectivities on its diagonal and F is the matrix of form factors. Thus, K looks like:

Using the form factor matrix computed before we can do matrix operations to compute 𝐾 which leads to this matrix:

 

The matrix K is diagonally dominant; thus, it is well suited to iterative methods as Gauss-Seidel, which is the method I used to solve the system of linear equations. 

The Gauss-Seidel Method is really fast with K matrix (only one iteration). Now to make it converge the spectral radius of K 𝜌(𝐾) must be less than 1. Therefore, the largest absolute value for an eigenvalue of K must be less than one.

Stimuli-Dependent 

There is a system of linear equations per stimuli required, so for the radiosity code I solved three systems of linear equations (red, green and blue).

Rendering

Once the radiosity values 𝐵 are known, I pass that information to vertex colors for each quadrilateral as shown in the next image:

 

The vertex values are used to bilinearly interpolate vertex colors when the rendering pass requests a color from the scene in a specific direction. 

These are vertex colors for face xy with z = 0 without interpolation. (I only render Upper Left Vertex Color for each Quad in this case): 

And this is face xz with y equal 0:

Now with interpolation respectively:

The radiosity solution is view-independent, therefore changing the view position does not need a new computation of radiosities nor form factors.

The Mesh is generated procedurally, for this image I used 10 elements per face side, each face side is 10 units long. 

Thank you!

If you want to look at the code:

Radiosity Code

Bibliography:

1. Hughes, J. F., & Dam, A. V. (2014). Computer Graphics Principles and Practice Third Edition (3rd ed.). Pearson Education.
2. Pharr, M., & Jakob, W. (2017). Physically Based Rendering From Theory to Implementation (3rd ed.). Cambridge, MA: Morgan Kaufmann.
3. Dutré, P., & Bala, K. (2006). Advanced Global Illumination (2nd. ed.). Wellesley, Massachusetts: A K Peters.
4. Neumann, L., Feda, M., Kopp, M., & Purgathofer, W. (1995). A New Stochastic Radiosity Method for Highly Complex Scenes. G. Sakas Et Al. (eds.), Photorealistic Rendering Techniques.
5. Cohen, M. F., & Wallace, J. R. (1993). Radiosity and Realistic Image Synthesis. Cambridge, MA: AP Professional.
6. Cohen, M. F., & Greenberg, D. P. (1985). The Hemi-Cube a radiosity solution for complex environments. Cornell University, Ithaca, N.Y.: ACM.
7. Strang, G. (2016). Introduction to Linear Algebra. (5th ed.). Wellesley, MA: WellesleyCambridge Press.
8. Stroustrup, B. (2013). The C++ Programming Language. (4th ed.). Pearson Education.