## Sampling microfacet BRDF

I’m working on microfacet brdf model for my renderer these days, noticing that it is more than necessary to provide a separate sampling method for microfacet brdf instead of using the default one, which is usually used for diffuse like surfaces and highly inefficient for brdf with spiky shape, such as mirror like surfaces. The following image is one generated by the default sampling method: The left monkey has pure reflection brdf which is mentioned in my previous blog post, the right one uses the microfacet model with zero as roughness value. I was expecting similar result for both of the monkeys while things turned out to be wrong here, we can barely see the reflection from the monkey. Actually there is nothing wrong, the fact is that the convergence rate of default sampling the microfacet model with 0 roughness is extremely low. As long as there are enough samples, it will reach the appearance similar to the left one. However the number of samples being enough can be arbitrarily high depending on how spiky your brdf is.

The right way of sampling microfacet brdf is mentioned in this paper. What I want to record in this blog is how these conclusions are derived from the original microfacet model. Using these better sampling method, we will get similar result for those two monkeys eventually.

### Why Default Sampling Method is Inefficient

So why is this default sampling method inefficient for mirror like brdf model. First of all, the default sampling method goes with the following pdf: $p_h(\omega) = \dfrac{cos(\theta)}{\pi}$

It achieves good result for diffuse surfaces like Lambert, OrenNayar. In fact, it is not a sampling method for diffuse brdfs, a sampling method for lambert brdf will respect a pdf with constant value, it also involves the cosine factor located in the rendering equation or LTE. ${L_{o}(\omega_{o}) = \int L_{i}(\omega_{i}) *f( \omega_{i}, \omega_{o} ) *cos(\theta_{i}) d\omega_{i}}$

Since it is extreme difficult, if not impossible, to have any knowledge on the incident radiance, the only thing we have is the brdf and cosine factor. The sampling method mentioned above could be efficient if the product of brdf and cosine factor is the dominant one in this equation. It works pretty well in practice most of the time.

However the efficiency of default sampling method drops quickly as the brdf becomes the dominant factor with roughness going to zero. For microfacet brdf here, its extreme form is like a dirac-delta function, it is almost impossible to get a sample with brdf having a non-zero value. And for the lucky ones that do have non-zero brdf values, it could reach super-high value due to its low probability, introducing high variance to the sampling result. In other words, the convergence rate is quite low, which is exactly why we are seeing the above image. It can be even treated as a bug. Obviously we need a better way to sample these microfacet brdfs.

### Microfacet Model

Microfacet Model is quite hot in real time rendering these years, it is the basics of physically based shading algorithm. Generally its basic form is like the following: $f(\omega_i,\omega_o,x) = \dfrac{F(\omega_i , h) G(\omega_i,\omega_o,h) D(h)}{4 cos(\theta_i) cos(\theta_o)}$

The first component is Fresnel, the second one is a shadowing factor called G term, the last one is normal distribution function (NDF). Comparing with the old ad-hoc bxdf model, microfacet model obeys the rule of energy conservation which allows the artist to change the roughness of the material through one parameter instead of two, specular color and specular power.

Among all those factors in this equation, the NDF term is usually the dominant one. The specific shape of NDF is largely affected by the roughness of the bxdf. In order to sample the microfacet brdf model, it is usually to sample NDF first getting a random microfacet normal that respects the NDF and then reflect the exitence radiance along the normal to generate the incident direction.

There are several kinds of NDF around. The following three will be covered in this blog and all that will be mentioned are isotropic ones.

• GGX
• Beckmann
• Blinn

Before we proceed, there is one rule that all NDF should follow: $\int_{\Omega} D(m) cos(\theta_m) d\omega = 1$

We will use this equation in our following derivation. Explaining the equation is outside the scope of this blog, however one can refer here for further detail.

#### Sampling GGX

Here is the basic form of GGX: $D(h) = \dfrac{\alpha^2}{\pi ((\alpha^2-1) cos^2\theta + 1 ) ^2}$

So the pdf respecting solid angle is like this: $p_h(\omega)=\dfrac{\alpha^2cos\theta}{\pi((\alpha^2-1) cos^2\theta+1)^2}$

Instead of sampling solid angle directly, we usually use spherical coordinate to sample it. So it is not the pdf respecting the solid angle that we are interested, it is the pdf respecting the spherical coordinates. $p_h(\theta,\phi)=\dfrac{\alpha^2cos\theta sin\theta}{\pi((\alpha^2-1) cos^2\theta+1)^2}$

The following task is quite straight forward, which is basically sampling according to a specific pdf. Inversion method goes pretty well for all NDFs in this blog. Notice that $\phi$ is not even in this equation, that said the NDF is totally isotropic and we can sample $\phi$ uniformly. For detail proof, please refer here. The only thing left here is how to sample $\theta$. In order to do that, we need to get the pdf for $\theta$ first. $p_h(\theta)=\int_{0}^{2\pi}p_h(\theta,\phi)d\phi = \dfrac{2\alpha^2cos\theta sin\theta}{((\alpha^2-1) cos^2\theta+1)^2}$

Let’s calculate the CDF next: $\begin{array} {lcl} P_h(\theta) & = & \int_{0}^{\theta} \dfrac{2\alpha^2 cos(t) sin(t)}{(cos^2t(\alpha^2-1)+1)^2}dt \\ & = & \int_{\theta}^{0} \dfrac{\alpha^2}{(cos^2t(\alpha^2-1)+1)^2}d(cos^2t) \\ & = & \dfrac{\alpha^2}{\alpha^2-1} \int_{0}^{\theta} d{\dfrac{1}{cos^2t(\alpha^2-1)+1}}\\ & = & \dfrac{\alpha^2}{\alpha^2-1} {(\dfrac{1}{cos^2\theta(\alpha^2-1)+1}-\dfrac{1}{\alpha^2})}\\ & = & \dfrac{\alpha^2}{cos^2\theta(\alpha^2-1)^2+(\alpha^2-1)}-\dfrac{1}{\alpha^2-1} \end{array}$

With a canonical random number $\epsilon$ equals to this CDF, we have the following equation: $\epsilon =\dfrac{\alpha^2}{cos^2\theta(\alpha^2-1)^2 + (\alpha^2-1)} - \dfrac{1}{\alpha^2-1}$

To solve this equation requires no more knowledge than mathematic in junior high school, I guess there is no need to show the whole process. The final solution of the above equation is: $\theta =arccos\sqrt{\dfrac{1-\epsilon}{\epsilon(\alpha^2-1)+1}}$ or $\theta =arctan(\alpha\sqrt{\dfrac{\epsilon}{1-\epsilon}})$

The above two equations are exactly the same thing. Choosing which one to use is just a matter of taste.

#### Sampling Beckmann

The derivation of sampling method for Beckmann is quite similar to the above procedure. Here is the Beckmann distribution: $D(h) =\dfrac{1}{\pi\alpha^2cos^4\theta}e^{-\dfrac{tan^2\theta}{\alpha^2}}$

The pdf respecting spherical coordinate is like this: $p_h(\theta,\phi)=\dfrac{sin\theta}{\pi\alpha^2cos^3\theta}e^{-\dfrac{tan^2\theta}{\alpha^2}}$

Again, $\phi$ can be sampled uniformly. The pdf for $\theta$ should be this: $p_h(\theta)=\dfrac{2sin\theta}{\alpha^2cos^3\theta}e^{-\dfrac{tan^2\theta}{\alpha^2}}$

The CDF for $\theta$ can be calculated this way: $\begin{array} {lcl} P_h(\theta) & = & \int_{0}^{\theta} \dfrac{2sin(t)}{\alpha^2cos^3t}e^{-\dfrac{tan^2t}{\alpha^2}}dt \\ & = & \int_{0}^{\theta} \dfrac{-2}{\alpha^2cos^3t}e^{-\dfrac{tan^2t}{\alpha^2}}d(cos(t))\\ & = &\int_{0}^{\theta} \dfrac{1}{\alpha^2}e^{-\dfrac{tan^2t}{\alpha^2}}d(\dfrac{1}{cos^2(t)}) \\ & = & \int_{0}^{\theta} \dfrac{1}{\alpha^2}e^{\dfrac{1}{\alpha^2}(1-\dfrac{1}{\cos^2t})}d(\dfrac{1}{cos^2(t)})\\ & = & \int_{\theta}^{0} d(e^{\dfrac{1}{\alpha^2}(1-\dfrac{1}{\cos^2t})})\\ & = & 1- e^{\dfrac{1}{\alpha^2}(1-\dfrac{1}{\cos^2\theta})} \end{array}$

Solving the equation of $P_h(\theta) = \epsilon$ gives the following solution: $\theta =arccos\sqrt{\dfrac{1}{1-\alpha^2 ln(1-\epsilon)}}$ or $\theta =arctan\sqrt{-\alpha^2\ln(1-\epsilon)}$

#### Sampling Blinn

Here is the Blinn NDF: $D(h)=\dfrac{\alpha+2}{2\pi}(cos\theta)^\alpha$

The pdf respecting to spherical coordinate is: $p_h(\theta,\phi)=\dfrac{\alpha+2}{2\pi}(cos\theta)^{\alpha+1}sin\theta$

Isolating $\phi$ gives the following: $p_h(\theta) = (\alpha+2)(cos\theta)^{\alpha+1}sin\theta$

This one is a lot simpler than the previous two, here is the CDF: $\begin{array} {lcl} P_h(\theta) & = & \int_{0}^{\theta} (\alpha+2)cos(t)^{\alpha+1}sin(t)d(t) \\ & = & \int_{\theta}^{0} (\alpha+2)cos(t)^{\alpha+1}d(cos(t)) \\ & = &\int_{\theta}^{0} d(cos(t)^{\alpha+2}) \\ & = &{1-cos(\theta)^{\alpha+2}} \end{array}$

Here is the relation between canonical random number and $\theta$ $\theta=arccos((1-\epsilon)^{\dfrac{1}{\alpha+2}})$

Since $\epsilon$ is a canonical random number, $1-\epsilon$ is also one. So we can simplify the above equation with the following one: $\theta=arccos(\epsilon^{\dfrac{1}{\alpha+2}})$

A little more about this sampling method. I’m not quite confident with this solution, although I can see nothing wrong in this derivation. PBRT gives a similar solution for Blinn which is: $\theta=arccos(\epsilon^{\dfrac{1}{\alpha+1}})$

And another open-source ray tracer, mitsuba, also adopts this way of sampling. I don’t quite understand the derivation in the book, so I’ll just stick to this one, which is also the one mentioned here. I tried the pbrt’s way of sampling, only minor difference can be noticed from the image.

#### One Extra Step

There is still one step further from being done. It is the incident direction sampling that we are interested, not the normal. The reason we are sampling the normal instead of incident direction directly is because NDF is usually the dominant factor in microfacet model. Generating the incident direction after sampling normal is a more efficient way than sampling the incident direction itself.

However the way we calculate PDF given an incident direction is different from the above ones respecting half-vector, whether it’s solid angle or spherical coordinate. What we have so far is the pdf for half-vector, a transformation is necessary. $p(\theta) = p_h(\theta) \dfrac{d\omega_h}{d\omega_i} = \dfrac{p_h(\theta)}{4 (\omega_o \cdot \omega_h)}$

### Conclusion

Here are the images generated after and before the new sampling method.  There are 64 samples per pixel and GGX is chosen as NDF, the three monkeys have different roughness values(0.0,0.5,1.0). As we can see from the image that the left-most monkey gets much better result than the default sampling method. And for the other two monkeys, we have similar result with the cosine-pdf sampling method.

There are already some research work improving the sampling method mentioned in this blog. I may need to try their method once I have some free cycles.

### Reference:

 Microfacet Models for Refraction through Rough Surfaces. https://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.html
 Dirac delta function. https://en.wikipedia.org/wiki/Dirac_delta_function
[
3] How is the NDF really defined. http://www.reedbeta.com/blog/2013/07/31/hows-the-ndf-really-defined/
[
4] Notes on importance sampling. http://blog.tobias-franke.eu/2014/03/30/notes_on_importance_sampling.html

### 29 thoughts on “Sampling microfacet BRDF”

1. Steve November 2, 2015 / 6:54 pm

You have a small typo in the equation at the end of the “Sampling GGX” section: the division by 1 – \epsilon should be inside the square root.

Liked by 1 person

• AGraphicsGuy November 3, 2015 / 1:02 am

Agreed, appreciated. It is already fixed.

Like

• Rafael June 28, 2020 / 8:12 am

How the “One Extra Step” formula is derived? How do you calculate dwh / dwi ?

Like

2. wallacekcy October 5, 2017 / 3:33 am

Hi,

Why does Beckmann’s pdf has a sine term instead of a cosine term like the others?

Thanks

Liked by 1 person

• AGraphicsGuy October 6, 2017 / 8:50 am

It should have a cosine term, I guess I missed it, :(. Already corrected.
Thanks for pointing it point.

Like

3. Alessandro Nardini November 3, 2017 / 3:12 pm

Very useful article, learning a lot from this! thank you for sharing!

Like

4. Fei Xia December 5, 2017 / 5:20 pm

it seems there exists small errors in the derivation of the cdf of both ggx and beckmann:
the second step of ggx and the fourth step of beckmann are missing the minus sign

Like

• Fei Xia December 5, 2017 / 5:26 pm

my bad, only the fourth step of beckmann is missing the minus sign.

Like

• AGraphicsGuy December 6, 2017 / 7:01 pm

I don’t see any issue of the derivation, I think there is no problem with the fourth step.

Like

5. Fei Xia December 5, 2017 / 5:41 pm

in my understanding the difference of result of bling is because you are sampling the normal, but the other is sampling directly the insident direction.
the answer live in this article:
https://developer.nvidia.com/gpugems/GPUGems3/gpugems3_ch20.html

Like

• AGraphicsGuy December 6, 2017 / 7:02 pm

I’m afraid that I don’t agree. It is normal that I sample in the three cases, there is no difference of bling comparing with others.

Like

6. g2-c02031016c3e5bf73dca2b8eb33f399a June 21, 2018 / 10:18 pm

Hi! Great post but I have a question: could you explain how you do this step in the calculation?

https://imgur.com/HJLWAFh

Look like magic…

Like

• AGraphicsGuy June 21, 2018 / 10:23 pm

Thanks.
I wish I can type latex here, which I don’t think so.
First of all, squared alpha in the nominator can be isolated because it is a constant. Try replace cos(t)^2 with x, you will get a new formula. And then what is left is very basic mathematics about integral calculus. I hope this helps.

Liked by 1 person

• g2-c02031016c3e5bf73dca2b8eb33f399a June 21, 2018 / 11:38 pm

Thanks. (but it’s very basic for you maybe, I still can’t properly integrate like you did, I’m heavily relying on Mathematica do to that for me 😀 )

Like

7. Aditya Sreekar August 18, 2018 / 12:48 pm

Awesome tutorial

Like

8. Hao Zhang March 8, 2019 / 2:59 pm

Hi! Great post, but I have a problem. In Sampling Part, for example, sampling GGX. ‘The PDF respecting solid angle’ you wrote is p(h,n,α)=D(h,n,α)|h⋅n|, and the solid angle means |h⋅n|, right? why should there be a |h⋅n|?

Like

• Hao Zhang March 8, 2019 / 3:03 pm

I figure out, after I read the link in the blog. Thanks.

Like

• AGraphicsGuy March 9, 2019 / 11:35 pm

To reply to your last question, none of the final conclusion reaches 0 when roughness goes to 0. Do notice that there is a curve from roughness to alpha, it is not a direct match. You can check out the code to get a better explanation of it.

Like

• Hao Zhang March 9, 2019 / 11:25 am

Thank you very much for your reply. I have one more question. When roughness equal to zero, the D(h) is always 0, doesn’t like a dirac-delta function which is 0 expect 1 point. So I can get good result at roughness close to 0, but it becomes totally black when the roughness is one. How did you deal with this 0 divided by 0 problem?

Like

9. Casper January 25, 2020 / 5:00 pm

Great blog, this is excactly the problem i’m encountering right now. I have one question though. For a faster converging image we need to sample in the direction of the reflection. In the end we have calculated a theta, that as far as I understand it, describes the angle between the normal and the outgoing ray. I’m not sure how this can be converted to an outgoing ray. Can we pick any outgoing direction which is at this angle?

Like

10. Hulda Erxleben October 15, 2020 / 5:22 pm

I’d like to thank you for the efforts you have put in penning this website. I am hoping to check out the same high-grade blog posts by you in the future as well. In truth, your creative writing abilities has motivated me to get my very own blog now 😉

Like

11. Jewel Heath October 17, 2020 / 4:44 am

Having read this I believed it was very enlightening. I appreciate you finding the time and energy to put this article together. I once again find myself spending a significant amount of time both reading and commenting. But so what, it was still worth it!

Like