Importance sampling microfacet GGX

Question

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))
}

Answer 1

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)
}