<p>Through help from several people, &amp; 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:</p>

<p><a href="https://i.sstatic.net/6mTVc.jpg" rel="noreferrer"><img src="https://i.sstatic.net/6mTVc.jpg" alt="enter image description here"></a></p>

<pre><code>// 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 &lt;= 0 || wi.Dot(wm) &lt;= 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)
}
</code></pre>

CC BY-SA 4.0
2020-08-25T00:36:12.330Z
Should `wo` in surface space or world space?
iaomw
CC BY-SA 4.0
2020-07-06T15:24:44.110Z
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)
vuoriov4
CC BY-SA 4.0
2018-05-23T14:09:10.840Z
2018-05-23T14:09:10.840Z
hunterloftis
<p>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 <a href="https://computergraphics.stackexchange.com/questions/130/trying-to-implement-microfacet-brdf-but-my-result-images-are-wrong">another user's question</a>, 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.</p>

<p>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).</p>

<p>I would really love to know how to translate <a href="https://hal.archives-ouvertes.fr/hal-01509746/document" rel="noreferrer">"A Simpler and Exact Sampling Routine for the GGX Distribution of Visible Normals"</a> into separate Sample() and PDF() functions.</p>

<p>Here's what I have for Lambert:</p>

<pre><code>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))
}
</code></pre>

<p>And here's my in-progress work for Microfacets:</p>

<pre><code>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)
}
</code></pre>

<p>The GGX BSDF relies on several physics functions I've also implemented:</p>

<pre><code>// 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))
}
</code></pre>

CC BY-SA 4.0
2018-05-22T08:24:14.637Z
`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
CC BY-SA 4.0
2018-05-22T08:44:52.740Z
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."
hunterloftis
CC BY-SA 4.0
2018-05-22T12:43:14.610Z
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
<raytracing><pathtracing>
Importance sampling microfacet GGX
hunterloftis
# Importance sampling microfacet GGX

(Asked by: username_0 on 2018-05-22)

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](https://computergraphics.stackexchange.com/questions/130/trying-to-implement-microfacet-brdf-but-my-result-images-are-wrong), 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"](https://hal.archives-ouvertes.fr/hal-01509746/document) 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))
}

```

**Tags:** raytracing, pathtracing

> **Comments:**
> **username_1**: `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  
> **username_0**: 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."  
> **username_1**: 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.  

---

## 1 answers

### Answer by username_0 on 2018-05-23. Score: 7

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:

[![enter image description here](https://i.sstatic.net/6mTVc.jpg)](https://i.sstatic.net/6mTVc.jpg)

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

```

> **Comments:**
> **username_3**: 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)  
> **username_2**: Should `wo` in surface space or world space?  

---

