# Importance sampling microfacet GGX

I have a renderer where the BxDF interface is Sample(), PDF(), and Eval(). The Lambertian BRDF is working well, and I believe I have properly implemented Eval for GGX based on another user's question, but I'm struggling to find a resource I can understand that explains how to generate samples, and then calculate the corresponding PDFs, for a GGX distribution.

I've found some of the original papers, but unlike Lambert the Sampling and PDF seem based on spherical coordinates (which I don't really understand - 3D Vectors are about the limit of my math).

I would really love to know how to translate "A Simpler and Exact Sampling Routine for the GGX Distribution of Visible Normals" into separate Sample() and PDF() functions.

Here's what I have for Lambert:

type Lambert struct {
R, G, B float64
}

func (l Lambert) Sample(out, normal geom.Direction, rnd *rand.Rand) geom.Direction {
return normal.RandHemiCos(rnd)
}

func (l Lambert) PDF(in, normal geom.Direction) float64 {
return in.Dot(normal) * math.Pi
}

func (l Lambert) Eval(in, out, normal geom.Direction) rgb.Energy {
return rgb.Energy{l.R, l.G, l.B}.Scaled(math.Pi * in.Dot(normal))
}


And here's my in-progress work for Microfacets:

func (m Microfacet) Sample(out, normal geom.Direction, rnd *rand.Rand) geom.Direction {
// TODO: better sampling
return normal.RandHemi(rnd)
}

func (m Microfacet) PDF(in, normal geom.Direction) float64 {
// TODO: PDF that matches a better sampling distribution
return 1 / (2 * math.Pi)
}

// https://computergraphics.stackexchange.com/questions/130/trying-to-implement-microfacet-brdf-but-my-result-images-are-wrong
// https://schuttejoe.github.io/post/ggximportancesamplingpart2/
func (m Microfacet) Eval(in, out, normal geom.Direction) rgb.Energy {
F := schlick2(in, normal, m.F0.Mean())  // The Fresnel function
D := ggx(in, out, normal, m.Roughness)  // The NDF (Normal Distribution Function)
G := smithGGX(out, normal, m.Roughness) // The Geometric Shadowing function
r := (F * D * G) / (4 * normal.Dot(in) * normal.Dot(out))
return m.F0.Scaled(r)
}


The GGX BSDF relies on several physics functions I've also implemented:

// Schlick's approximation of Fresnel
func schlick2(in, normal geom.Direction, f0 float64) float64 {
return f0 + (1-f0)*math.Pow(1-normal.Dot(in), 5)
}

// GGX Normal Distribution Function
// http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html
func ggx(in, out, normal geom.Direction, roughness float64) float64 {
m := in.Half(out)
a := roughness * roughness
nm2 := math.Pow(normal.Dot(m), 2)
return (a * a) / (math.Pi * math.Pow(nm2*(a*a-1)+1, 2))
}

// Smith geometric shadowing for a GGX distribution
// http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html
func smithGGX(out, normal geom.Direction, roughness float64) float64 {
a := roughness * roughness
nv := normal.Dot(out)
return (2 * nv) / (nv + math.Sqrt(a*a+(1-a*a)*nv*nv))
}

• based on spherical coordinates (which I don't really understand - 3D Vectors are about the limit of my math) you should really look into this - once you get the hang of it, it seems pretty obvious. with cartesian coordinates, you simply describe the extend on each axis for a point. with spherical coordinates, you rather specify the direction, in which the point lies. imagine you look straight up. the polar angle determines, how far you rotate your head down (y axis, y is up). the azimuthal angle determines, how far you rotate your body around (around the y-axis). that's the gist of it
– Tare
Commented May 22, 2018 at 8:24
• I understand the gist of it, but I don't understand it well enough to translate the referenced paper into an implementation. For example, I don't expect to understand this anytime soon: "We sample the projected area by generating samples on the half disks proportionally to their respective projected areas. For this, we use a polar parameterization of the disk (r, φ) in which we scale the angle φ in order to account for the difference of projected areas of its two halves." Commented May 22, 2018 at 8:44
• While I can't explain this to you atm, they are probably working with solid angles. This, too, is a concept you should familiarize yourself with when working with BxDFs. It is somewhat confusing, especially in the beginning, but it is really worth it. Basically you cut out a piece of a surface of a sphere and see how big the area of the cut out piece is in relation to the squared radius of the sphere. I doubt you can get to the bottom of your question without understanding solid angles first.
– Tare
Commented May 22, 2018 at 12:43

## 1 Answer

Through help from several people, & referencing and re-referencing the commented urls, I managed to get a sampling scheme with matching PDF that handles everything from perfect mirror surfaces to roughness = 1. Here's the eventual result:

// Cook-Torrance microfacet model
type Microfacet struct {
F0        rgb.Energy
Roughness float64
}

// https://schuttejoe.github.io/post/ggximportancesamplingpart1/
// https://agraphicsguy.wordpress.com/2015/11/01/sampling-microfacet-brdf/
func (m Microfacet) Sample(wo geom.Direction, rnd *rand.Rand) geom.Direction {
r0 := rnd.Float64()
r1 := rnd.Float64()
a := m.Roughness * m.Roughness
a2 := a * a
theta := math.Acos(math.Sqrt((1 - r0) / ((a2-1)*r0 + 1)))
phi := 2 * math.Pi * r1
x := math.Sin(theta) * math.Cos(phi)
y := math.Cos(theta)
z := math.Sin(theta) * math.Sin(phi)
wm := geom.Vector3{x, y, z}.Unit()
wi := wo.Reflect2(wm)
return wi
}

// https://schuttejoe.github.io/post/ggximportancesamplingpart1/
// https://agraphicsguy.wordpress.com/2015/11/01/sampling-microfacet-brdf/
// https://en.wikipedia.org/wiki/List_of_common_coordinate_transformations#From_Cartesian_coordinates_2
func (m Microfacet) PDF(wi, wo geom.Direction) float64 {
wg := geom.Up
wm := wo.Half(wi)
a := m.Roughness * m.Roughness
a2 := a * a
cosTheta := wg.Dot(wm)
exp := (a2-1)*cosTheta*cosTheta + 1
D := a2 / (math.Pi * exp * exp)
return (D * wm.Dot(wg)) / (4 * wo.Dot(wm))
}

// http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.html
func (m Microfacet) Eval(wi, wo geom.Direction) rgb.Energy {
wg := geom.Up
wm := wo.Half(wi)
if wi.Y <= 0 || wi.Dot(wm) <= 0 {
return rgb.Energy{0, 0, 0}
}
F := fresnelSchlick(wi, wg, m.F0.Mean()) // The Fresnel function
D := ggx(wi, wo, wg, m.Roughness)        // The NDF (Normal Distribution Function)
G := smithGGX(wo, wg, m.Roughness)       // The Geometric Shadowing function
r := (F * D * G) / (4 * wg.Dot(wi) * wg.Dot(wo))
return m.F0.Scaled(r)
}

• For anyone using this, i would like to point out that (x, y, z) are not absolute coordinates. I got the correct result by generating an orthonormal basis where y = surface normal. (This is mentioned in schutte's article though) Commented Jul 6, 2020 at 15:24
• Should wo in surface space or world space? Commented Aug 25, 2020 at 0:36