Pre-filtered environment map, deriving the equation

Question

I'm reading through this article, and more specifically I'm trying to derive the equation that would explain the implementation the following shader (still in the same article):

#version 330 core
out vec4 FragColor;
in vec3 localPos;

uniform samplerCube environmentMap;
uniform float roughness;

const float PI = 3.14159265359;

float RadicalInverse_VdC(uint bits);
vec2 Hammersley(uint i, uint N);
vec3 ImportanceSampleGGX(vec2 Xi, vec3 N, float roughness);

void main()
{       
    vec3 N = normalize(localPos);    
    vec3 R = N;
    vec3 V = R;

    const uint SAMPLE_COUNT = 1024u;
    float totalWeight = 0.0;   
    vec3 prefilteredColor = vec3(0.0);     
    for(uint i = 0u; i < SAMPLE_COUNT; ++i)
    {
        vec2 Xi = Hammersley(i, SAMPLE_COUNT);
        vec3 H  = ImportanceSampleGGX(Xi, N, roughness);
        vec3 L  = normalize(2.0 * dot(V, H) * H - V);

        float NdotL = max(dot(N, L), 0.0);
        if(NdotL > 0.0)
        {
            prefilteredColor += texture(environmentMap, L).rgb * NdotL;
            totalWeight      += NdotL;
        }
    }
    prefilteredColor = prefilteredColor / totalWeight;

    FragColor = vec4(prefilteredColor, 1.0);
}  

Which is my understanding suppose to implement the computation of the integral

$$ I = \int_{\Omega} L(\omega_i)d\omega_i $$

However I'm not entirely sure since NdotL is also used in the computation so I might be wrong.

It seems to me the idea is given $N$ to sample a random halfway vector $H$, and based on this we can construct, by reflection, the light direction vector $L$, which would be random as well by construction. The article also mentions that $H$ is sampled from a normal distribution.

But anyway suppose I have samples $H^{1} ,\ldots, H^{N}$ and therefore $L^{1},\ldots, L^{N}$ samples, the latter I suppose would define $\omega_{i}^{1} \ldots, \omega_{i}^{N}$ (solid angles).

By importance sampling

$$ I \approx \frac{1}{N} \sum_{j=1}^{N} \frac{L(\omega_i^{j})}{pdf(\omega_i^{j})} \;\;\;\; (1) $$

If what I wrote is correct than I don't know where the dot product and the normalization factor

$$ W = \sum_{j=1}^{N} N\cdot L^{j} $$

comes from, namely how do I go from $(1)$ to

$$ I \approx \frac{\sum_{j=1}^{N}L(\omega_i^{j}) N \cdot L^{j}} {\sum_{j=1}^{N} N \cdot L^{j}} $$

Update: I'm pretty sure many of you know this already, but the shader above is for a pre-filtered environment map. Which means in this context the actual integral should be given by

$$ I(N) = \int_{\Omega} L(\omega_i) N \cdot L d\omega $$

This explains the factor $N \cdot L$ in the shader, the distribution function used is the GGX distribution, I assume plugging this stuff together might give me the formula.

Just updating because it might give a better clue.