Microfacet BRDFs from Scratch

In this article, I will derive these famous Specular and Diffuse BRDFs from scratch.

BRDFspec=F(ωi,ωh) G1(ωi,ωo,ωh) D(ωh)4 ∣ωg⋅ωi∣ ∣ωg⋅ωo∣\boxed{\text{BRDF}_{\text{spec}} = \frac{F(\omega_i, \omega_h), G_1(\omega_i, \omega_o, \omega_h), D(\omega_h)}{4,|\omega_g \cdot \omega_i|,|\omega_g \cdot \omega_o|}}BRDFspec​=4∣ωg​⋅ωi​∣∣ωg​⋅ωo​∣F(ωi​,ωh​)G1​(ωi​,ωo​,ωh​)D(ωh​)​​

BRDFdiff=1∣ωg⋅ωi∣ ∣ωg⋅ωo∣∫Ω1π ⟨ωi⋅ωm⟩ G1(ωi,ωo,ωm) D(ωm) ⟨ωo⋅ωm⟩ dωm\text{BRDF}{\text{diff}} = \frac{1}{|\omega_g \cdot \omega_i|,|\omega_g \cdot \omega_o|} \int\Omega \frac{1}{\pi}, \langle \omega_i \cdot \omega_m \rangle, G_1(\omega_i, \omega_o, \omega_m), D(\omega_m), \langle \omega_o \cdot \omega_m \rangle, d\omega_mBRDFdiff​=∣ωg​⋅ωi​∣∣ωg​⋅ωo​∣1​∫Ω​π1​⟨ωi​⋅ωm​⟩G1​(ωi​,ωo​,ωm​)D(ωm​)⟨ωo​⋅ωm​⟩dωm​

I will try to make the derivations as intuitive as possible with figures and explanations. The interactive HTML sections have been written by AI (The free version of Claude) as I don't know how to do frontend HTML programming. The other figures were made by me using inkscape.

This derivation is not novel, I read through 5 very important papers shuffling all of them at once, and I got the idea while going through them, to create a straight-forward article for someone like me, who is interested in mathematics and derivations of equations, but for whom, certain jumps in derivations don't come so naturally.

​dωd\omegadωDifferential solid angleD(ωm)D(\omega_m)D(ωm​)Normal Distribution Function (NDF)G1G_1G1​Masking-shadowing functionFFFFresnel factor⟨ωi⋅ωm⟩\langle \omega_i \cdot \omega_m \rangle⟨ωi​⋅ωm​⟩Clamped dot product. The BRDF​

When a light hits any "thing", it will either be transmitted (refracted), reflected, or absorbed. Absorption happens when the atom that is hit by the light takes the energy of the light and it reaches an excited state. Eventually, the energy is spit out, if it's emitted as visible light then that phenomenon is known as Photoluminescence. This isn't really important to Rendering as we rarely have to simulate that phenomena on production or real-time renderers. BRDF deals with the reflection of light. If we want to emulate the transmission as well, it'll require another function with a very similar name, the BTDF.

Bi-directional Reflectance Distribution Function (BRDF) is a scary name, but it is quite simple. It gets a bad reputation because when people try to learn it, they get met with very scary looking equations (some of which we will derive in this article).

The BRDF is simply the ratio of outgoing radiance in a given direction to the differential irradiance from a specific incoming direction. Radiance is the power of light in a point in a single direction. Irradiance is the power of light in a point from all directions. That is why we use differential irradiance. If we take an integral of the differential irradiance, we get the whole irradiance. I came across this cheatsheet in the graphics programming discord for important terms in rendering such as radiance, irradiance and what not.

So, BRDF =Lo/Ei= L_o / E_i=Lo​/Ei​ then? Yes! But, we don't actually know the EiE_iEi​ in rendering. We will have to figure out the Irradiance by taking the incoming light from each angle around the point and integrating it. Taking the irradiance of each angle is not practical at all, there's infinite angles as there are infinite numbers. So we use something called a differential solid angle, which is an infinitesmally small surface 'patch' in a sphere. As that patch is extremely small, the irradiance along the surface at any point will be virtually the same, so we can measure the irradiance at one point in the surface and then project that patch to the point.

The surface area of the patch, despite being extremely small, is larger than that of a singular 'point' of whose BRDF we are measuring. Depending on the angle at which the patch is located with respect to the normal, the incoming light gets spread out and misses the point.

Figure 1. Interactive visualiser showing how the same solid angle dωd\omegadω projects onto a larger surface element dAdAdA as the incident angle θ\thetaθ increases.

Due to this, the incoming light at a point through a solid angle will then become:

Licos⁡θi dωiL_i \cos\theta_i , d\omega_iLi​cosθi​dωi​

LiL_iLi​ is the incoming light from a direction iii, θi\theta_iθi​ is the light incident angle wrt normal, and dωd\omegadω is the differentiable solid angle.

BRDF=LoLicos⁡θi dωi\text{BRDF} = \frac{L_o}{L_i \cos\theta_i , d\omega_i}BRDF=Li​cosθi​dωi​Lo​​

In rendering, we mostly want to figure out the Radiance to a specific directions, so we write it as:

Conservation of energy — the total reflected energy cannot exceed the incident energy:

∫Ωfrcos⁡θo dωo≤1\int_\Omega f_r \cos\theta_o , d\omega_o \leq 1∫Ω​fr​cosθo​dωo​≤1

Reciprocity — swapping incoming and outgoing directions leaves the BRDF unchanged:

fr(θi,ϕi,θo,ϕo)=fr(θo,ϕo,θi,ϕi)f_r(\theta_i, \phi_i, \theta_o, \phi_o) = f_r(\theta_o, \phi_o, \theta_i, \phi_i)fr​(θi​,ϕi​,θo​,ϕo​)=fr​(θo​,ϕo​,θi​,ϕi​)

A 'smooth' surface to our eye will have extreme difference in elevation as we zoom in. There will be many bumps in the microscopic scale.

Figure 2. A macrosurface that appears smooth has microsurfaces with varying normals at the microscopic scale.

If we want a realistic simulation of light then we cannot ignore this scale. We are not going to nanogeometry however, we will only think of any surface as being built up of very small microsurfaces. Each point we want to sample on the macrosurface (which may appear as a pixel on our screen) will have microsurfaces. So if we want to figure out the total radiated light in a certain direction for a certain point in the macrosurface, we need to find the radiation for that direction throughout all the points in the microsurface and sum all of it. Since these microsurfaces are very small, the out-radiance angle for each point of the microsurface will be parallel to one-another. This can be visualised as each point projecting itself to a straight line.

Figure 3. The microsurface points all project onto a line of width cos⁡θωo\cos\theta_{\omega_o}cosθωo​​ in the outgoing direction ωo\omega_oωo​.

The points in the microsurface will have different normals so each point's out-radiance will be spread out in the outgoing direction. If the normal of a point in the microsurface is parallel to the outgoing direction, the out-radiance will be projected 1:1 to the projected 'line', while for any other angle the out-radiance will be spread out.

This formula acts as a weighted average of the out-radiance of each point of the microsurface. An intuitive metaphor is to imagine a folding phone, if it's folded 90 degrees and you're looking at it. The folded screen that is directly facing your eyes will contribute more to the brightness than the one that is facing up, and is 90 degrees to your eyes.

Figure 4. The folding phone analogy: the screen directly facing the viewer contributes more brightness than the one at 90 degrees. Each microsurface point contributes proportionally to its projected area.

For a macrosurface with microsurfaces, let us create a function D(ω)D(\omega)D(ω) which counts the amount of points in the microsurface with a specific normal ω\omegaω.

δω\delta_\omegaδω​ is known as a dirac delta, think of it as a function that discards what is not expected. δω\delta_\omegaδω​ means that the delta will discard every input which isn't ω\omegaω. ωm(pm)\omega_m(p_m)ωm​(pm​) gives the normal of the microsurface point pmp_mpm​. So, D(ω)D(\omega)D(ω) gives the total points in the microsurface with the given normal ω\omegaω.

There's an infinite amount of vectors as normals that a point can have. These normals exist in a unit sphere as well. A unit sphere also has all the normals in 3-D. For each point in the microsurface with a normal ωm\omega_mωm​, there's a point in the unit sphere with that normal. For a subset of point in the microsurface M′⊆M\mathcal{M}' \subseteq \mathcal{M}M′⊆M, there is a subset of unit sphere Ω′⊆Ω\Omega' \subseteq \OmegaΩ′⊆Ω with the same normals.

Integrating for each point in M′\mathcal{M}'M′ gives the same 'area' as integrating for each point in a unit sphere with the same normals, because D(ωm)D(\omega_m)D(ωm​) gives the weight for each of those normals. Reminder that even in the unit sphere integral, the D(ωm)D(\omega_m)D(ωm​) goes through each point in the microsurface and checks the normal of each point. So, for the whole microsurface, or the whole unit circle, it gives the microsurface area.

Since a surface is made up of microsurfaces (Figure 2), the projections will be that of each point in the microsurfaces. Not all the points in the microsurfaces are actually visible, depending on the viewing angle. The visibility also depends on what's known as a 'microsurface profile'. The same distribution of normals can have different behaviors and have different visibility depending on how the normals are distributed along the microsurfaces.

Figure 5. Two microsurfaces with the same normal distribution but different profiles can have different visibility. Some points are occluded depending on the viewing direction ωo\omega_oωo​.

The projected area of the microsurfaces can be viewed as such:

GGG here is a function that evaluates to 0 if a point pmp_mpm​ is blocked, and 1 if it's not blocked. Clamp function is used to cull out backfacing microgeometries, as it'll clamp negative value to 0, so if there's a point which faces the opposite of ωo\omega_oωo​, it'll not contribute to the projected area. We can create another GGG function that takes ωm\omega_mωm​ instead of pmp_mpm​ [2]. For that:

Here, D(ωm)D(\omega_m)D(ωm​) gives weight to the specific ωm\omega_mωm​. Meaning, if there's more points with a specific ωm\omega_mωm​ which are also visible, their contribution has more weight. Since the projected area is the area of all visible points in the microsurface, for a surface of unit area, we know that its area from (ii). From (ii) and (v) we get:

In equation (iii), let us replace the denominator with (ii). We can replace the whole integral in the numerator with another numerator that also takes the area (v), and for the function L(ωo,pm)L(\omega_o, p_m)L(ωo​,pm​), which is a function of the point, we can just convert it to a function of ωm\omega_mωm​: L(ωo,ωm)L(\omega_o, \omega_m)L(ωo​,ωm​). (The reason it's so easy to convert is because the radiance depends on the normal of the point regardless).

Here, Dωo(ωm)D_{\omega_o}(\omega_m)Dωo​​(ωm​) acts as a distribution of the visible normals, it gives value between 0 and 1.

L(ωo,ωm)L(\omega_o, \omega_m)L(ωo​,ωm​) is the outgoing radiance at a direction ωm\omega_mωm​. Using basic calculus, we can rewrite it as:

Basically, we differentiate the function wrt dωid\omega_idωi​ and integrate it along dωid\omega_idωi​. From (i), we can say that the micro BRDF, BRDFm\text{BRDF}_mBRDFm​ (It is a function of ωo\omega_oωo​, ωi\omega_iωi​, and ωm\omega_mωm​) for a specific microsurface is:

We can ignore the back faces using the clamp, so:

Now, let's differentiate equation (viii) with respect to the incoming irradiance:

We replace dLodL_odLo​ from (xi) in (xiv):

The macro BRDF is given in (i). We can replace the dLodL_odLo​ with (xv), then replace Dωo(ωm)D_{\omega_o}(\omega_m)Dωo​​(ωm​) from (ix). cos⁡(θωo)\cos(\theta_{\omega_o})cos(θωo​​) and cos⁡(θi)\cos(\theta_i)cos(θi​) can be taken out of the integral and written as ∣ωg⋅ωi∣|\omega_g \cdot \omega_i|∣ωg​⋅ωi​∣ and ∣ωg⋅ωo∣|\omega_g \cdot \omega_o|∣ωg​⋅ωo​∣ as they depend on the macro surface normal:

This statistical model of BRDF covers the light before it bounces from the microsurface. After it bounces, it could get blocked by another microsurface, for that, we use a funct