A Little bit of Background

When we are doing physically-based rendering, a direct path-tracing approach is sampling many points from the light source and extending a light ray starting from each point. Rays interact (refraction or reflection) with the objects in the scene and some eventually reach the camera. We sum the final radiance of rays within each pixel of the camera and obtain the resulting image.

In order to simulate how a material affects the radiance coming towards it, we have the Bidirectional Reflectance Distribution Function (BRDF) model. It tells us how is an incidence light ray reflected on the surface of a material. In order to keep consistent with the paper (Jakob et al. 2014 ) which introduces the discrete microfacet theory used for this project, I will use the definition from the paper:

Plainly speaking, the result of rendering is just the total amount of light projected onto each pixel of our camera. In other words, we want to obtain the power scattered on a finite area around the incidence point towards a finite solid angle around the outgoing direction, and we integrate the BRDF defined above like so:

Another important point to consider is the scale at which we evaluate our BRDF. We can treat each geometry in the scene as a whole (the macro scale) or we can take a closer look at a material to account for the characteristics it has at this level (the meso scale). However, for example, if we are considering a material like the cover of a MacBook, both of these scales will produce results that say it is perfectly smooth while in reality, we know that this is not true. Therefore, we have to see the material through a microscope (the micro scale) to evaluate its properties correctly. This is where the microfacet theory comes into place.

The Microfacet Theory

At the micro scale, we assume that the surface consists of many facets (represented mathematically by normals pointing in different directions). A light ray would only contribute to the final radiance if it was reflected to the path hitting the camera given the facets. In other words, the direction of the reflected ray is dependent on the directions of the facets. Therefore, we need to find out how many of the facets are pointing in the “right“ directions. It is certainly unrealistic to model the facets one by one, as it would be incredibly inefficient. However, we can make use of statistical physics and say our facets follow the Beckmann distribution:

We can use other distributions such as the Blinn distribution or the GGX distribution as well, but in this project, I only used the Beckmann distribution.

The microfacet model, as a BRDF, acts as the foundation for rendering glints. Its equation is:

The Fresnel reflection coefficient F tells us how much light is reflected and transmitted on an optical object. The shadowing and masking factor G takes into account that a ray which is reflected after hitting a facet might get occluded by another facet. G is still a statistical approximation which is specific to the microfacet distribution used. For Beckmann distribution, we can express G as:

The D term is continuous; hence the materials rendered will look very smooth and uniform. However, as we can see earlier, glints appear to be more “grainy“, which is nothing like the example given below. Is this still the same thing? How can we achieve the glinty effect instead? The trick here is to replace the continuous D with a discrete one which we will elaborate in the next section.

The Discrete Microfacet Distribution

As you might have guessed, the glinty effect uses a discrete D term instead of a continuous one:

If we substitute the discrete D term into the microfacet BRDF, substitute the resulting BRDF to the integral and evaluate the function with approximation, we get the final version of the power function we want to implement in this project:

where D hat is a multiscale version of the microfacet distribution. Since this blog emphasises the implementation of the theory, details of how the final version is derived are omitted. If you are interested, please refer to section 3.1 of (Jakob et al. 2014).

The key point is to find an algorithm that can evaluate D hat efficiently.

Though the original paper (Jakob et al. 2014) has introduced two algorithms of evaluation, they are both too complicated for the scope of this project. Therefore, an improved algorithm (Atanasov and Koylazov, 2016) is used instead.

There are two major steps in the algorithm:

Build and traverse a quad-tree for a given incidence ray;
Generate the exact UV coordinates for the leaf nodes of the quadtree and evaluate D hat within the footprint parallelogram A.

Step I: Build and Traverse a Quad-Tree

When we export the UV texture map of a shape from Blender, the map should be a unit square with bottom-left as (0,0) and top-right as (1,1). We first define the number of flakes N we want to sprinkle on this map, and we draw a footprint parallelogram A on it with ray differentials. Then we split the square (i.e. node) into 4 equal-sized squares and distribute N to each of the sub-square with a multinomial distribution of probabilities p=(1/4, 1/4, 1/4, 1/4) (note that the distribution is done on each split, and we need to make sure that the sum of the flakes in sub-squares is equal to the number of flakes in their parent). If a sub-square contains no flake or does not overlap A, it is immediately discarded. We keep splitting a node until either it contains less than K flakes or the stack depth T is reached. The D hat calculation is done in the leaf nodes.

This “split and visit” pattern is combining quad-tree building and traversal. We use depth-first search to traverse the tree.

Efficient Approximation of Multinomial Distribution

Although the number of flakes in each node should be drawn according to the multinomial distribution, the brute force implementation will be slow during runtime because it is an expensive calculation. Therefore, we only compute an approximation to the original distribution using the multidimensional normal distribution for a node if it has more than n flakes (in Jakob et al. 2014, n=8). The multidimensional normal distribution has a mean of μ and covariance Σ.

Mathematically, given two uniform random samples, a pair of independent random variables with a standard normal distribution (mean 0, variance 1) can be generated by Box-Muller transform:

We need two pairs of such samples and put them in a matrix Z where:

We want to transform this matrix to have a mean of μ:

and a 4×4 covariance matrix Σ:

We can get the resulting vector by:

Note that the samples drawn from the distribution are not usually integers, so we need a way to convert them to the right type. More details will be given in the implementation section.

Step II: Evaluate D Hat

While traversing through the quad-tree, we also evaluate D hat if the node is a leaf. We first randomly generate the exact UV coordinates for each flake inside the node. If 𝐴 contains a flake's coordinates, we generate a random normal 𝑚_𝑗 for that flake according to the desired microfacet distribution (in our case it's the Beckmann distribution). Then we reflect our incidence ray 𝜔_𝑖 at 𝑚_𝑗 to get a per-flake outgoing direction 𝑟_𝑗. 𝑟_𝑗 associates a specific cone with radius 𝛾 (in radians). If the cone contains the reflected ray ω_o of ω_i, we compute the flake weight 𝑤_𝑗=𝑚_𝑗⋅𝑟_𝑗 which approximates the microfacet distribution. Then D hat is simply:

Implementation

Framework Nori

Before I dive into my implementation, I’ll shortly introduce the framework upon which my project is built.

Nori is an open-source educational rendering engine developed by Dr Wenzel Jakob in 2012. It includes a basic C++ framework upon which my project is constructed. More specifically, my job is extending the features of the framework, such as adding more BRDF models and types of integrators.

If you are interested in more details, please click on the link provided in the title.

Overview

The glinty effect is only meaningful under direct illumination with a small light source because otherwise, it behaves like a normal microfacet BRDF. Therefore, two integrators are extended which are whitted (whitted.cpp) and multiple-importance sampling path tracing (path_mis.cpp). The former is mainly used for testing. The integrators were written in previous projects and they are not the focus of this project; hence they will only be briefly explained in later sections. A class Glints (in glints.cpp) is created to extend BRDF. The only difference between Glints and Microfacet is that Glints has this extra function called evalGlint which makes use of the microfacet BRDF with the D hat. The diffuse term in evalGlint is removed because it is not meaningful to the glinty effect.

The actual D hat computation happens in class QuadTree (in quadtree.h). The description of each function inside this class is given below:

evalD(): Given a footprint, the UV map keeps getting split into sub-squares and the sub-squares that overlap with the footprint is pushed onto a stack. If a node has larger than n flakes, it is split; otherwise, it is used for D hat evaluation.
getBboxForFootprint(): Since we need to check if a node overlaps with the footprint, we do this by computing the bounding box of it and see if the bounding box overlaps with that of the node. Note that this might be inefficient if the footprint is a highly skewed parallelogram.
distributeFlakesMultinomial(): Generate four samples from a multinomial distribution with the Box-Muller transform.
distributeFlakesNaive(): Generate four samples from a multinomial distribution with brute force.

For more details on the implementation, please check out the source code at the bottom of this blog.

How to Get the Pixel Footprint in UV Space

There are many ways to define a tiny area in the UV space as the unit of the D hat evaluation. The method used in the paper is ray differentials. To understand this, we need to consider the three spaces for rendering: the pixel space, the world space, and the UV (i.e. texture) space. The question is to figure out how is the area of one pixel projected onto the UV space. Both Physically Based Rendering: From Theory to Implementation (PBRT) and Mitsuba2 have full implementation of ray differentials and hence the code is used for this project. (Mitsuba3 was still in the pre-release version at the time I was doing this project. Nonetheless, the implementation for ray differentials is the same as that of Mitsuba2)

A point within a pixel is sampled in the renderBlock() function in the main.cpp file. The traceRayDifferentials() is called on the sample to get the origins and directions of the two rays offset by 1 pixel to the right and up separately. Then when we call scene->rayIntersect() inside an integrator given a ray and an empty intersection object its, the function computes the values for dp_du and dp_dv as attributes of its. Later on, if we encounter an object with glinty BRDF, we call its.compute_uv_partials() to get duv_dx and duv_dy. Now we know how the UV coordinates change with respect to the pixel coordinates, the footprint is simply:

// The footprint vertices are in clockswise order with top-left being the first one.
Point2f footprint[4] = {};
footprint[0] = its.uv + its.duv_dx * -.5f + its.duv_dy * .5f;
footprint[1] = its.uv + its.duv_dx * .5f + its.duv_dy * .5f;
footprint[2] = its.uv + its.duv_dx * .5f + its.duv_dy * -.5f;
footprint[3] = its.uv + its.duv_dx * -.5f + its.duv_dy * -.5f;

How to Distribute the Flakes to each Sub-Square

The brute force approach (distributeFlakesNaive()) is straightforward. We initialise a vector representing four bins and loop through the number of elements n we want to generate. A bin index is generated with a uniform random number generator and the selected bin gets +1 sample.

The Box-Muller transform (distributeFlakesMultinomial()) is less intuitive. To get the non-integer vector result, we implement the equation given previously. Then we cast the vector as int and randomly increase or decrease the number of samples in a bin until the sum of all bins is equal to n.

How to Approximate D Hat

When we get a leaf node, we loop through the current number of flakes inside it and generate random UV coordinates for a flake. If a flake is inside the footprint, we generate a random normal 𝑚_𝑗 for that flake according to Beckmann distribution (Warp::squareToBeckmann()) and reflect the incident ray 𝜔_𝑖 at 𝑚_𝑗 to get a per-flake outgoing direction 𝑟_𝑗. If the cosine angle between 𝜔_𝑜 and 𝑟_𝑗 is greater than or equal cos(𝛾), we add the flake weight 𝑤_𝑗=𝑚_𝑗⋅𝑟_𝑗 to the accumulative float D. Note that 𝑚_𝑗⋅𝑟_𝑗 is only an approximation of the local microfacet distribution and can be replaced by Warp::squareToBeckmannPdf() to have more accurate results. Eventually, D hat =𝐷/𝑁 is returned.

Since we need millions of flakes to see the glinty effect and it takes way too long to render on my local machine, I set a WEIGHT_SCALE parameter to increase the contribution of each valid flake so that I do not need to render that many flakes to see a result.

How to Ensure D Hat is Consistent

We are modelling a discrete microfacet distribution and technically the locations of the flakes and their normal vectors should be consistent all the time. The magic of this algorithm is that we do not need to store all that information to evaluate D hat, we just need to make sure that the random number generator generates the same set of numbers for a node if needed.

Each node of the quadtree is given an ID with the type of size_t. The ID for the root node is 1 and the ID for a sub-node is computed as the ID of its parent bit shifted to the left by 2 units plus the position index of that node (the order of children goes in the top left, top right, bottom left, and bottom right). For example, the top right child of the root node has an ID of (1 << 2) + 1 = 101 (the ID is binary). In this case, each child will have a unique and consistent ID that is used to seed the random number generator and get the same sequence of numbers every time. Since all calculations that need random samples use the same numbers, the results are consistent. In this case, we ensure that the underlying discrete microfacet distribution remains unchanged.

Integrators

Integrators trace the path travelled by a single input light ray a limited number of times with Russian roulette and output the final radiance of it. They are called " integrators" because theoretically ray tracing is math integrations. Although a few integrators were implemented in past projects, only whitted.cpp and Multiple Importance Sampling path_mis.cpp are used. whitted.cpp is a simplified version of an integrator and only takes direct illumination into account (as you can see from the Cornell Box scenes in the validation section, the ceilings are black because they are not lit up), which is ideal for test cases. On the other hand, path_mis.cpp considers direct and indirect illumination. The idea is that every time the ray intersects in the scene we produce two light paths:

A path directly connects this intersection to the emitter (if the path is not blocked by an obstacle);
A sampled outgoing path that does not hit an emitter (so that we do not double count the contribution of the emitter).

Note that when the ray hits a specular material, its outgoing direction is fixed. Therefore, we do not connect this type of intersection by emitter sampling which considers more than one point on the emitters. Instead, we account for the contribution of the emitter at this intersection with path tracing only. For diffuse materials, we weigh and combine the contribution of both the emitter sampling and BRDF sampling.

Below is the pseudo-code for the MIS path tracer:

if (there is no emitter in the current scene) {
    return Color3f(0.0f);
}

initialise all necessary variables;

for (int i = 0; i < maxDepth; ++i) {
    if (there is no intersection for the current ray) {
        return L;
    }

    // specularBounce returns true if the last intersection is specular
    if (the current intersection is on an emitter) {
        if (i == 0 || specularBounce) {
            add the contribution of the current intersection directly to L;
        }
        else {
            add the weighted contribution of the current intersection to L;
        }
    }

    if (the current intersection is on a diffuse material) {
        sample a point from emitters;
        if (the current material requires ray differentials) {
            // Microfacet is a type of diffuse material that needs ray differentials
            compute and store the ray differentials in its fields;
        }
        add the weighted contribution of the emitter point to L;
    }

    sample an outgoing direction from the current intersection;

    update the variables; 

    // Russian roulette termination criteria
    if (i > 3) {
        // force ray to bounce at least 3 times
        calculate continuation probability;
        if (the uniform random variable > continuation probability) {
            break;
        }
        else {
            throughput /= continuation probability;
        }
    }
}

return L;

This integrator is used for rendering the final image for the project.

Validation

Unfortunately, there is no implementation of glints in either PBRT or Mitsuba2. Extra tests need to be conducted to validate the implementation.

The Ray Differentials are Correct

Even though the code for ray differentials is taken from Mitsuba2, it is still necessary to check whether it is doing the right thing. Correct ray differentials do not change the distribution of the render result as the conversion produces correct UV coordinates centred around the intersection. However, small differences between different unwrapping schemes are accepted. As mentioned by Atanasov and Koylazov, 2016, oftentimes texture coordinates have bad quality; hence we cannot expect different unwrapping to give us the exact same result.

The section below shows the comparison of the results between Unwrapping I and Unwrapping II with all other parameters (i.e. the roughness parameter alpha and the number of flakes) being the same.

As we can see, Unwrapping II is much more “dense” than Unwrapping I so the flakes are more centred at the peak of the distribution. Nonetheless, two unwrapping methods produce the same distribution. Therefore, the ray differentials are doing the right thing.

D is Consistent

The best way to test whether the underlying D hat is consistent is to see whether rendering twice the same scene (with the exact same parameters) produces the exact same image.

Below is the light source at the same position (and angle) rendered twice. As you can see the two resulting images are identical.

Another way to test whether the D hat is consistent is to see whether the most flake-concentrated area (i.e. the distribution itself) follows the light source. Below are the images rendered with the light source at different angles. Indeed, the flake distribution is coherent with the light source.

Glints Converge to Continuous Microfacet

As the number of flakes on the object increases, the glint BRDF should get closer and closer to a continuous microfacet distribution.

As we can see the more flakes there are, the more glint BRDF looks like a continuous microfacet distribution. Therefore the implementation is correct.

My Final Submission

Recall that the final project has a theme called “loop“. Well, now let’s see what I have done.

*The list of BRDFs I used:*
*Background: Diffuse Texture Mapping*
*Hearts, Igea (the statue in the centre), Bimba (two statues of kids on the side): Glint BRDF*
*Lucy (statues surrounding Igea), Balls, Roses: Continuous Microfacet*

I take my inspiration from the Kaleidoscope since it repeatedly loops through beautiful patterns as you turn. The most prominent feature of the Kaleidoscope is that the patterns are symmetrical. I then build the scene in Blender with some classic 3D models like so:

Note that all objects are inside a Cornell box but you cannot see it from the final picture because the walls only act as emitters (i.e. light sources) and reflective surfaces to light up the scene.

Next, I want the colour scheme to match the sparkling girly style I envision when implementing the glints; hence I choose the rosy vintage colour palettes which are included in the final.xml file.

View the Source Code

Below is the portal to my GitHub gist.

Click to Open

Acknowledgements

My biggest challenge in this project was that I was the only student in the class implementing the glint effect. Furthermore, the TAs could not help me because this was an advanced topic not in their research areas. Therefore, I had to turn up for Dr Wenzel Jakob’s limited office hours every Thursday with the list of questions prepared in advance and prioritised. Dr Jakob was extremely patient and clearly explained everything to me. I highly appreciate his help. Without it, I wouldn’t be able to overcome this challenge.

References

Matt Pharr, Wenzel Jakob, and Greg Humphreys. 2016. Physically Based Rendering: From Theory to Implementation (3rd ed.). Morgan Kaufmann Publishers Inc., San Francisco, CA, USA. https://www.pbr-book.org/
Jakob, W., Haˇsan, M., Yan, L.-Q., Lawrence, J., Ramamoorthi, R., and Marschner, S. (2014). Discrete stochastic microfacet models. ACM Trans. Graph., 33(4).
Atanasov, A. and Koylazov, V. (2016). A practical stochastic algorithm for rendering mirror-like flakes. In ACM SIGGRAPH 2016 Talks, SIGGRAPH ’16, New York, NY, USA. Association for Computing Machinery.
Merlin Nimier-David, Delio Vicini, Tizian Zeltner, and Wenzel Jakob. 2019. Mitsuba 2: a retargetable forward and inverse renderer. ACM Trans. Graph. 38, 6, Article 203 (December 2019), 17 pages. https://doi.org/10.1145/3355089.3356498
Wenzel Jakob, Sébastien Speierer, Nicolas Roussel, and Delio Vicini. 2022. DR.JIT: a just-in-time compiler for differentiable rendering. ACM Trans. Graph. 41, 4, Article 124 (July 2022), 19 pages. https://doi.org/10.1145/3528223.3530099

C++ Photo Realistic Rendering Glint Effect: Discrete Microfacet with Stochastic Algorithm

What are glints?

A Little bit of Background

The Microfacet Theory