Importance sampling GGX NDF - fireflies and bright final result https://schuttejoe.github.io/post/ggximportancesamplingpart1/ - I have problem implementing this method. Somebody asked question about this before ( Can&#39;t understand the Importance sampling GGX ) but it doesn't fully answer my question/problem.
Above blog post describe importance sampling method based on importance sampling NDF. Let's start with Cook-Torrance BRDF:
$$
f_r = \frac{F(w_i, w_m) ~ G_2(w_i, w_o, w_m) ~ D(w_m)}{4 ~ |w_i \cdot w_g| ~ |w_o \cdot w_g|}
$$
After finding PDF for Normal Distribution Function and applying
$$
|w_i \cdot w_g|
$$
term from rendering equation + found PDF, we have final equation:
$$
\frac{F(w_i, w_m) ~ G_2(w_i, w_o, w_m) ~ |w_o \cdot w_m|}{|w_o \cdot w_g| ~ |w_m \cdot w_g|}
$$
I'm adding final equation just for completion sake. All transformations are rather long and it's better to check blog post that I've sent. At the end of post, there is code that should be working correctly. However, it doesn't give correct results in my case and my suspicion is that I'm filling gaps in a wrong way:
I'm using Correlated Multi-Jittered Sampling here:
float SmithGGXMaskingShadowing(float3 n, float3 l, float3 v, float a2)
{
    float dotNL = saturate(dot(n, l));
    float dotNV = saturate(dot(n, v));

    float denomA = dotNV * sqrt(a2 + (1.0f - a2) * dotNL * dotNL);
    float denomB = dotNL * sqrt(a2 + (1.0f - a2) * dotNV * dotNV);

    return 2.0f * dotNL * dotNV / (denomA + denomB);
}

float a2 = roughness * roughness;
float theta = acos(sqrt((1.0f - brdfSample.x) / ((a2 - 1.0f) * brdfSample.x + 1.0f)));
float phi = 2.0f * PI * brdfSample.y;

float3 normalTS = float3(0, 1, 0);

float3 wo = -incomingRayDirTS;
float3 wm = float3(sin(theta) * cos(phi), cos(theta), sin(theta) * sin(phi));
float3 wi = 2.0f * dot(wo, wm) * wm - wo;

float3 F = Specular_F_Schlick(specularAlbedo.rgb, saturate(dot(wi, wm)));
float G = SmithGGXMaskingShadowing(wm, wi, wo, a2);
float weight = abs(dot(wo, wm) / (dot(normalTS, wo) * dot(normalTS, wm)));

rayDirWS = normalize(mul(wi, tangentToWorld));

if (dot(normalTS, wi) &gt; 0.0f &amp;&amp; dot(wi, wm) &gt; 0.0f)
    throughput = F * G * weight;
else
    throughput = 0.0f;

My main concern is - am I using normal in tangent space (w_g) in a correct way? Schutte Joe mentions in a comment that:

// -- Ensure our sample is in the upper hemisphere
// -- Since we are in tangent space with a y-up coordinate
// -- system BsdfNDot(wi) simply returns wi.y

That's why I decided to use w_g (or normalTS) as float3(0, 1, 0).
Here is how it looks for 1024 samples with path length equal 4:

In comparison, here is method based on visible normals described here - https://schuttejoe.github.io/post/ggximportancesamplingpart2/ ; For 1024 samples, path length 4, results are much better. However, based on data provide by author of the post, there shouldn't be that much difference:


Edit: @B_old is right, I mixed TS and WS. I present corrected results above -   image is slightly better but still has very high variance and some bright spots.
I've made a few changes in my code, so let's start with basic image that I had problem we that we'll state as a problem at the beginning:

After changing part of the code from:
float weight = abs(dot(wo, wm) / (dot(normalTS, wo) * dot(normalTS, wm)));
if (dot(normalTS, wi) &gt; 0.0f &amp;&amp; dot(wi, wm) &gt; 0.0f)

To this (i.e. moving part of calculations to world-space):
float weight = abs(dot(wo, wm) / (dot(triangleNormal, incomingRayDirWS) * dot(normalTS, wm)));
if (dot(triangleNormal, rayDirWS) &gt; 0.0f &amp;&amp; dot(wi, wm) &gt; 0.0f)

I am getting better results. Here are series of images with different settings:
1024 spp, 4 bounces

5000 spp, 4 bounces

5000 spp, 1 bounce

Notice that 4 bounces version is having more fireflies and noise. Long specular paths might be a reason of introducing noise. Also, 5000 spp are getting darker than 1000 spp. Reason is that path is probably introducing no changes, therefore average color is getting darker. You can see it in the last line of main post. I'm going to put it below for reference:
if (dot(triangleNormal, rayDirWS) &gt; 0.0f &amp;&amp; dot(wi, wm) &gt; 0.0f)
    throughput *= F * G * weight;
else
    throughput *= 0;

Last thing. Here is image generated with better method ( https://schuttejoe.github.io/post/ggximportancesamplingpart2/ ). 1024 spp, 4 bounces:


Overall, I got rid of weird noise and artifacts. Algorithm got closer to better method and I think it can be treated as fixed.