Questions about implementing glTF 2.0 material system in path tracer

[summer coral]

I am trying to implement the [glTF 2.0 material system][1] in my (monte carlo) path tracer. I am stuck on a few details, and need assistance:

0. In the specification, the vectors L and V are used which typically refer to the vector from the intersection point towards a light source and the vector from the intersection point towards the camera respectively.

I do not have direct illumination implemented yet, thus there are no dedicated light sources in my scenes. The only "light source" is a background color which could be seen as "ambient lighting". Thus, I have chosen to interpret L as wi in my implementation where wi is the sampled bounce ray. My motivation for this is that I am bouncing in that direction in order to find out what the light contribution from that angle is. Thus, in some way L and wi are the same thing(ish). I do not know if this is a correct thing to do.

Furthermore, since there are many rays bouncing around there isn't a "given viewing direction"; rather, I am choosing to interpret V as wo where wo is the incoming ray direction. Again, I do not know if this is a correct thing to do.

1. Figuring out how to sample a new bounce direction in SpecularBRDF::sample, as well as what the PDF is.

The specification uses a microfacet BRDF that uses some complicated functions. Evaluating the BRDF is easy enough; however, it is not trivial to figure out how to sample a new direction as well as what the PDF is. At least not for me.

2. Figuring out how to sample a new bounce direction in Dielectric::sample, as well as what the PDF is.

Similar problem as with 2. above. A Dielectric material has a SpecularBRDF and a DiffuseBRDF. I have no clue how I am supposed to take into account both of them when wanting to sample a new direction as well as defining the PDF.

Below follows my work in progress(!!!) implementation of the material system mentioned above.

[sharp aspen]

Your understanding of L and V as the incoming ray and sample directions is correct.

BRDF sampling is quite involved. Perhaps these will help:
https://www.youtube.com/watch?v=NA-sKoByYlI
https://pbr-book.org/4ed/Reflection_Models/Roughness_Using_Microfacet_Theory
https://agraphicsguynotes.com/posts/sample_microfacet_brdf/
https://schuttejoe.github.io/post/ggximportancesamplingpart1/
https://jcgt.org/published/0007/04/01/paper.pdf

[onyx jay]

look at vndf sampling for importance sampling the microfacet normal, then you can just do reflect(wo, microfacet_normal) to get wi, the pdf is ggx_D()*G_1/(4.0*dot(wo, geometric_normal)

for dielectrics you should have 2 sampling strategies, one for diffuse and one for specular then choose between the 2 using the fresnel term as the probability of sampling the specular strategy (it doesnt have to be the fresnel term but it does a little bit of importance sampling)

the pdf of that is just specular_pdf*F + (1.0 - F)*diffuse_pdf

(just like how the bsdf should be evaluated)

[onyx jay]

also i found the slides explaining it to be way easier to understand than the paper so here it is https://jcgt.org/published/0007/04/01/slides.pdf

[summer coral]

I have some questions.

So, first and foremost. During sampling, what exactly is alpha_x and alpha_y? I only have alpha. My intuition is telling me that it has to do with anisotropy, which in my case would mean that alpha_x = alpha_y = alpha since I am not taking into account anisotropic materials (for now at least).

Given that, the sampling process is simply copying the code-snippet from the slides. That's straight forward enough.

---

I am a bit confused regarding the calculation of the PDF. In the attached image you see the microfacet BRDF specified in the glTF 2.0 spec. In the other image you see the code-snippet for sampling included in the slides. Note the comment above the function mentioning the PDF.

1. I have a feeling that the Smith joint masking-shadowingfunction G(form glTF 2.0 spec) is NOT equal to the G1 function which you include in the PDF. Is that correct? What exactly is G1? Why 1, is there a G2?

2. Is the microfacet distribution function D from the glTF 2.0 spec the same as ggx_Dthat you mention?

3. What are the input arguments to ggx_D and G1 in the case of evaluating the PDF? E.g. the glTF 2.0 spec uses N and H where N = geometric_normal and H = normalize(wi + wo). However, I have a feeling that those are not what are to be input to ggx_D in this case? Regarding G1, I don't even know what function that is so I can't really say much about what the inputs might be. I am imagining something like this...

const glm::dvec3 H = glm::normalize(wi + wo);  // wi = glm::reflect(wo, Ne) where Ne = microfacet normal
const double ggx_D = D(Ne, H);  // Ne = microfacet normal
const double G1 = G_1(Ne, wo, H);
const double pdf = ggx_D * G1 / (4.0 * glm::max(glm::dot(wo, normal), 0.0001));

I don't know if this is correct.

4. You said PDF = ggx_D() * G_1 / ( 4.0 * dot(wo, geometric_normal); however, The paper states PDF D_Ve(Ne) = G1(Ve) * max(0, dot(Ve, Ne)) * D(Ne) / Ve.z.

Here is what I have so far:

[[nodiscard]] Sample sample(const glm::dvec3& wo, const glm::dvec3& normal) const override {
        // source: https://jcgt.org/published/0007/04/01/slides.pdf
        
        const double U1 = Util::RandomDouble();
        const double U2 = Util::RandomDouble();
        const double alpha_x = alpha;
        const double alpha_y = alpha;

        // transforming view direction to the hemisphere configuration
        glm::dvec3 Vh = glm::normalize(glm::dvec3((alpha_x * wo.x, alpha_y * wo.y, wo.z)));

        // creating an orthonormal basis
        const glm::dvec3 T1 = (Vh.z < 0.9999) ? glm::normalize(glm::cross(glm::dvec3(0.0, 0.0, 1.0), Vh)) : glm::dvec3(1.0, 0.0, 0.0);
        const glm::dvec3 T2 = glm::cross(Vh, T1);

        // parametrization of the projected area
        const double r = glm::sqrt(U1);
        const double phi = 2.0 * Util::PI * U2;
        const double t1 = r * glm::cos(phi);
        double t2 = r * glm::sin(phi);
        const double s = 0.5 * (1.0 + Vh.z);
        t2 = (1.0 - s)*glm::sqrt(1.0 - t1*t1) + s*t2;

        // reprojection onto hemisphere
        const glm::dvec3 Nh = t1*T1 + t2*T2 + glm::sqrt(glm::max(0.0, 1.0 - t1*t1 - t2*t2))*Vh;

        // transforming the normal back to the ellipsoid configuration (to microfacet normal)
        const glm::dvec3 Ne = glm::normalize(glm::dvec3(alpha_x * Nh.x, alpha_y * Nh.y, glm::max(0.0, Nh.z)));

        // sample new direction
        const glm::dvec3 wi = glm::reflect(-wo, Ne);

        // PDF calculation
        // TODO: ...
        const double pdf = 0.0;

        return {wi, pdf};
    }

... and here is the code for V, G and D

    [[nodiscard]] double Chi(const double x) const {
        if (std::fabs(x) <= 0.0) {
            return 1.0;
        }
        return 0.0;
    };

    // Smith joint masking-shadowing function
    [[nodiscard]] double G(const glm::dvec3& wi, const glm::dvec3& wo, const glm::dvec3& N) const {

        const glm::dvec3 H = glm::normalize(wi + wo);

        const double NdotWI = glm::max(glm::dot(N, wi), 0.0001);
        const double NdotWO = glm::max(glm::dot(N, wo), 0.0001);

        const double HdotWI = glm::dot(H, wi);
        const double HdotWO = glm::dot(H, wo);

        const double numeratorL = 2.0 * NdotWI * Chi(HdotWI);
        const double denomL = NdotWI + glm::sqrt(alpha2 + (1 - alpha2) * NdotWI * NdotWI);
        const double quotientL = numeratorL / denomL;

        const double numeratorR = 2.0 * NdotWO * Chi(HdotWO);
        const double denomR = NdotWO + glm::sqrt(alpha2 + (1 - alpha2) * NdotWO * NdotWO);
        const double quotientR = numeratorR / denomR;

        return quotientL * quotientR;
        
    }

    // GGX - Normal Distribution Function
    [[nodiscard]] double D(const glm::dvec3 N, const glm::dvec3 H) const {
        const double NdotH = glm::max(glm::dot(N, H), 0.0001);

        const double numerator = alpha2 * Chi(NdotH);
        const double temp = (NdotH * NdotH * (alpha2 - 1) + 1);
        const double denom = Util::PI * temp * temp;

        return numerator / denom;
    }

    [[nodiscard]] double V(const glm::dvec3& wi, const glm::dvec3& wo, const glm::dvec3& N) const {

        const double numerator = G(wi, wo, N);

        const double NdotWI = glm::max(glm::dot(N, wi), 0.0001);
        const double NdotWO = glm::max(glm::dot(N, wo), 0.0001);
        const double denom = 4.0 * NdotWI * NdotWO;

        return numerator / denom;
    }

[summer coral]

I've been trying to figure stuff out and I believe that I have made some progress related to point 1) from above. I found these https://ubm-twvideo01.s3.amazonaws.com/o1/vault/gdc2017/Presentations/Hammon_Earl_PBR_Diffuse_Lighting.pdf slides from Earl Hammon which shed some light on some stuff. The math is really hairy though.

The "Smith joint masking-shadowing" function G from the glTF 2.0 specification seems to take on something of the formG2(wi,wo, *) = G1(wi,*) * G1(wo,*) where G1 in this case takes the form of the expression seen in the image below. That expression is specifically for G1(wi,*).

Thus, the answer to the point 1) should be that the G_1 that you are referring to is the expression in the image below. However, it is still unclear what the arguments to it should be. I have a feeling it is the calculated microfacet normal Ne (which we reflect around to produce our sampled direction wi), and... well I don't know what the second argument should be?

The number n in Gn seems to refer to the number of directions we're taking into account regarding geometric visibility. Using G2 = G1*G1.

I don't know if this understanding regarding point 1) is correct.

something like this instead of the above version of G.

    [[nodiscard]] double G_1(const glm::dvec3& N, const glm::dvec3& A, const glm::dvec3& H) const {
        const double NdotA = glm::max(glm::dot(N, A), 0.0001);
        const double HdotA = glm::dot(H, A);

        const double numerator = 2.0 * NdotA * Chi(HdotA);
        const double denom = NdotA + glm::sqrt(alpha2 + (1 - alpha2) * NdotA * NdotA);

        return numerator / denom;
    }

    [[nodiscard]] double G(const glm::dvec3& wi, const glm::dvec3& wo, const glm::dvec3& N) const {

        const glm::dvec3 H = glm::normalize(wi + wo);

        const double NdotWI = glm::max(glm::dot(N, wi), 0.0001);
        const double NdotWO = glm::max(glm::dot(N, wo), 0.0001);

        const double HdotWI = glm::dot(H, wi);
        const double HdotWO = glm::dot(H, wo);

        const double G1_wi = G_1(N, wi, H);
        const double G1_wo = G_1(N, wo, H);

        return G1_wi * G1_wo;
    }

[summer coral]

I am a bit confused on how to do this. In the glTF 2.0 spec, the Fresnel term is calculated as seen in the image below. Notice how it uses VdotH where V=wo in my case and H = normalize(wi + wo). However, how am I supposed to evaluate this given that I haven't yet sampled a bounce direction wi? That's the whole point, that I want to use the Fresnel term as a statistical random variable to help me make the decision about whether to sample diffuse or specular. 🤔

... and here is the incomplete WIP material system https://pastebin.com/XNRfwLvC.

The following scene was rendered with only indirect lighting (no direct illumination nor emissive objects). The supposedly Metal material is coming out all black 😄 obviously SpecularBRDF class is not correct, likely the SpecularBRDF::sample(...) function is wrong. Everything else in the scene is a Dielectric; however, it technically is not correct either because I am not probabilistically sampling the DiffuseBRDF or SpecularBRDF. It is just hard-coded to always sample the DiffuseBRDF so I guess it can just be seen as a pure Lambertian material for now. At least that looks good and correct 😄 but again, technically incorrect because the Dielectric material should take into account the SpecularBRDF too 🙂

Ignore the blue square at the top, it'll eventually become an area light or something like that 😄

[onyx jay]

Yeah in isotropic materials alpha x == alpha y.

1. yes G1 isn't the G term, the G term is G1(wo)*G1(wi), so G1 is the masking function (tells you the probability a ray doesn't intersect a microfacet in some direction. I'm on mobile so I can't post the code for it but it's here https://schuttejoe.github.io/post/ggximportancesamplingpart2/

2. I haven't seen the gLTF spec so idk for sure but it probably is

3. the inputs of D term is the same as the one used when evaluating the bsdf, G1 has the view vector as the input

4. that's just the pdf of sampling the normal iirc, you also need to divide by the jacobian of reflecting the ray

[onyx jay]

Well the fresnel when calculating the probability can just be an approximate so I would just use the geometric normal for calculating fresnel when picking between strategies and use microfacet normal when actually sampling the bsdf, but there might be a better method

[summer coral]

Yeah, I've been trying to figure things out but I can't seem to get stuff to work properly. The current state of the material system is this https://pastebin.com/YEFQvZvM. The result is seemingly too bright/washed and overly grainy/noisy. The goal is to implement the glTF 2.0 mat system, which seems quite straight forward... but is proving to be problematic (or I just suck lol). https://registry.khronos.org/glTF/specs/2.0/glTF-2.0.html#appendix-b-brdf-implementation

@sharp aspen @onyx jay my bad if it's not appropriate to ping like this, do you have any ideas what could be wrong? see pastebin link above for the material system code. It would also be interesting to see other people's implementations to compare with, you don't happen to have your own path tracer? 😄

[sharp aspen]

Not of the kind you're after. The thing about learning these things is that nobody will ever read your code to help you out. It's just too big of an ask. At best people will engage with a couple of lines of code if there is a very specific and isolated question you have.

The best way for you to get what you want is to really understand the basics of the model you are implementing. And that can be done by independent research of the literature. You may get lucky and find a very similar implementation on github but it will be readable only when you know the theory of sampling, mixing etc. I don't think there is a good way around reading up on these things for a couple of weeks.

Also, it's better to focus on one aspect of your problem at a time. For example, instead of trying to research how all of the gltf model works, try to get just brdf sampling first. It also makes for smaller questions and more likely help on discord

[summer coral]

fair enough, I completely get it

[sharp aspen]

if (metalness == 0.0f) {
  //--------------------- Diffuse ------------------------

  vec2 Xi = getRandomVec2(rngState);

  sampleDir = importanceSampleCosine(Xi, N);

  /*
      The discrete Riemann sum for the lighting equation is
      1/N * Σ(brdf(l, v) * L(l) * dot(l, n)) / pdf(l))
      Lambertian BRDF is c/PI and the pdf for cosine sampling is dot(l, n)/PI
      PI term and dot products cancel out leaving just c * L(l)
  */
  localCol = albedo;

} else {
  //--------------------- Specular ------------------------

  vec2 Xi = getRandomVec2(rngState);
  // Get a random halfway vector around the surface normal (in world space)
  vec3 H = importanceSampleGGX(Xi, N, roughness);

  // Generate sample direction as view ray reflected around h (note sign)
  sampleDir = normalize(reflect(-V, H));

  float NdotL = dot_c(N, sampleDir);
  float NdotV = dot_c(N, V);
  float NdotH = dot_c(N, H);
  float VdotH = dot_c(V, H);

  vec3 F = fresnel(VdotH, F0);
  float G = smiths(NdotV, NdotL, roughness);

  /*

      The following can be simplified as the D term and many dot products cancel out

      float D = distribution(NdotH, roughness);

      // Cook-Torrance BRDF
      vec3 brdfS = D * F * G / max(0.0001, (4.0 * NdotV * NdotL));

      float pdfSpecular = (D * NdotH) / (4.0 * VdotH);
      vec3 specular = (L(sampleDir) * brdfS * NdotL) / pdfSpecular;

  */

  // Simplified from the above

  localCol = (F * G * VdotH) / (NdotV * NdotH);
}
col *= localCol + emissive;