Ray Tracing, why does the GGX make my rendering result darker and darker for each frame? I am trying to add GGX and PBR textures to my GPU ray tracer. The result is very strange, it's getting darker for each frame.


        float ao = texAO.sample(textureSampler, hitRecord.uv).r;
        
        float3 albedo = texAlbedo.sample(textureSampler, hitRecord.uv).rgb;
        albedo = pow(albedo, 2.2); //albedo = float3(1.0);
        
        float metallic = texMetallic.sample(textureSampler, hitRecord.uv).r;
        float3 tNormal = texNormal.sample(textureSampler, hitRecord.uv).xyz;
        tNormal = normalize(tNormal * 2 - 1);
        
        float roughness = texRoughness.sample(textureSampler, hitRecord.uv).r;
        roughness = max(roughness, 0.001);
        
        float3 nx, ny;
        CoordinateSystem(normal, nx, ny);
        float3x3 stw = { nx, ny, normal };
        
        normal = stw * tNormal;
        CoordinateSystem(normal, nx, ny);
        stw = { nx, ny, normal };
        
        float3 wi, wh;
        float3 wo = transpose(stw) * (ray.direction);
    
        auto pSpecular = 1/(2-metallic);
        
        auto r0 = randomF(seed);
        auto r1 = randomF(seed);
        
        if (randomF(seed) &gt; pSpecular) {
            
            float sinTheta = sqrt(r0);
            float cosTheta = sqrt(1-r0);
            
            float phi = 2.0 * M_PI_F * r1;
            
            float x = sinTheta * cos(phi);
            float y = sinTheta * sin(phi);
            float z = cosTheta;
            
            wi = { x, y, z };
            wh = normalize(wo + wi); // not specular
            
        } else {
        
            auto a1 = roughness;
            auto a2 = a1 * a1;
            
            auto theta = acos(sqrt((1-r0) / ((a2-1)*r0+1)));
            auto phi = 2 * M_PI_F * r1;
            
            auto x = sin(theta) * cos(phi);
            auto y = sin(theta) * sin(phi);
            auto z = cos(theta) * 1;
            
            //wo = transpose(stw) * (ray.direction);
            wh = normalize( {x, y, z} );
            wi = reflect(wo, wh);
        }
        
        if (wi.z &lt;= 0 || dot(wh, wi) &lt;= 0) { return false; }
    
        auto D = NDF(wh, roughness);
        auto G = GX(wo, wi, roughness);
        auto F = FX(wi, wh, albedo, metallic);
        
        float denominator = 4.0f * abs(wo.z * wi.z);
        
        auto specular = D * G * F / max(FLT_MIN, denominator);
        
        auto diffuse = albedo / M_PI_F;
        
        auto kS = F;
        auto kD = (1 - metallic) * (1.0f - kS);

        ray = Ray(hitRecord.p, stw * wi);
        scatRecord.attenuation = (kD * diffuse + specular) * ao;
        
        //auto pdf = PDF(wh, roughness);
        //scatRecord.attenuation /= abs(pdf);
        
        return true;

// Normal Distribution Function
float NDF(const thread float3&amp; h, float roughness) {
  //  GGX / Trowbridge-Reitz
  
  float a1 = roughness;
  float a2 = a1 * a1;
  float NoH = abs(h.z);
  
  float d = (NoH * a2 - NoH) * NoH + 1;
  return a2 / (d * d * M_PI_F);
}

// Geometry Function
float GX(const thread float3&amp; wo, const thread float3&amp; wi, float roughness) {
  // Schlick, remap roughness and k

  if (wo.z &gt;= 0 || wi.z &lt;= 0)
      return 0;

  float k = pow(roughness + 1, 2) / 8.f;
  
  k = pow(roughness, 4);
  
  float G1_wo = wo.z / (wo.z*(1 - k) + k);
  float G1_wi = wi.z / (wi.z*(1 - k) + k);
  return G1_wo * G1_wi;
}

// Fresnel
float3 FX(const thread float3&amp; wi, const thread float3&amp; h, const thread float3 albedo, float metallic) {
  // Schlick’s approximation
  // use a Spherical Gaussian approximation to replace the power.
  // slightly more efficient to calculate and the difference is imperceptible
  
  float3 F0 = mix(float3(0.04f), albedo, metallic);
  float HoWi = dot(h, wi);
  return F0 + (float3(1.0f) - F0) * pow(2.0f, (-5.55473f * HoWi - 6.98316f) * HoWi);
}
http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.htm - here is reference with different geometry shading equations
http://jcgt.org/published/0003/02/03/paper.pdf - paper explaining how G term works and on what assumptions it is build.

Looking again at your code, I'm not sure (based on first link to Karis blog) if your G term is calculated correctly.
  float k = pow(roughness + 1, 2) / 8.f; 
  k = pow(roughness, 4);

You're overwritting k with roughness^4, which doesn't fall into any Schlick equation that can be seen in Karis blog and if I recall correctly, I haven't seen that in any papers. Another note is, that you're passing ωo and ωi to your geometry function, but based on your comment that this is Schlick G and overall context, I think that you should use NoL and NoV (dot first pointToLight direction and then viewDir).
I was trying to do my own image, but I found [Blinn 1977] paper which sums up how Smith G model works (take a look at figure 1) - https://www.microsoft.com/en-us/research/wp-content/uploads/1977/01/p192-blinn.pdf#page=2&amp;zoom=auto,-264,459
Some time ago, I was playing a little with different D, F, G equations - https://github.com/komilll/dxFramework/blob/master/dxFramework/Shaders/PS_BRDF.hlsl - take a look at PS_BRDF.hlsl and PS_BRDF_Helper.hlsl for further reference. I was basing my equations on Karis blog, some are modified base on other papers e.g. Frostbite and Disney diffuse are based on respective papers if I recall correctly.

I'm not sure if that's causing a problem, but it's worth fixing for sure. Also I don't think that you're using for summing frames
(previous_average * frame_count + current_color) / (frame_count + 1) 

is correct, neither for uniform or importance GGX sampling. You should derive pdf in that case, it can be found in link that I've sent you before - http://cwyman.org/code/dxrTutors/tutors/Tutor14/tutorial14.md.html - it's generally nice tutorial, with good step by step introduction to DXR.