I don't understand why the PBR implementation of Metorpolis Light Transport scales the final image by b / mutationsPerPixel.

The authors write:

Each Metropolis iteration within "Run nChains Markov chains in parallel" has splatted contributions with weighted unit luminance to the Film so that the final average film luminance before Film::WriteImage() is exactly equal to mutationsPerPixel.

(see http://www.pbr-book.org/3ed-2018/Light_Transport_III_Bidirectional_Methods/Metropolis_Light_Transport.html at the very end of the page.)

For simplicity, assume that we are considering a fixed path length $k\in\mathbb N$ and a single fixed $(s,t)$-strategy. Formally, let

• $M$ denote the scene surface set, $\mathcal B(M)$ denote the Borel $\sigma$-algebra on $M$ and $\sigma_M$ denote the surface measure on $\mathcal B(M)$
• $E:=M^{\{0,\:\ldots\:,\:,k\}}$, $\mathcal E:={\mathcal B(M)}^{\otimes\{0,\:\ldots\:,\:,k\}}$ and $\lambda:=\sigma_M^{\otimes\{0,\:\ldots\:,\:,k\}}$
• $f:E\to[0,\infty)^3$ denote the measurement contribution function
• $q$ be the probability density on $(E,\mathcal E,\lambda)$ corresponding to the $(s,t)$-strategy satisfying $\{q=0\}\subseteq\{f=0\}$

Once more, for simplicity, assume that the generated $(E,\mathcal E)$-valued chain $(X_n)_{n\in\mathbb N}$ is independent and $$X_n\sim q\lambda\;\;\;\text{for all }n\in\mathbb N.$$ If it's easier to understand for you, you might want to take a look at the importance_sampling_integrator::render() method in this code, where I've implemented these simplifications.

Now let's consider the measurement of a single pixel value $$I:=\int hf\:{\rm d}\lambda,$$ where $h:E\to[0,\infty)$ is the image reconstruction filter corresponding to this pixel. Assume (as PBRT does) that $h$ is a Box filter with radius 1/2. So, if $\psi:E\to R$ is the canonical mapping from the path space $E$ to the raster space $R=[0,a)\times[0,b)$, where $a\in\mathbb N$ and $b\in\mathbb N$ are the vertical and horizontal resolution of the image, respectively, and the pixel we've fixed is pixel $(i,j)\in J:=\{0,\ldots,a-1\}\times\{0,\ldots,b-1\}$ then $h$ is simply the indicator function of the set $$B:=\left\{x\in E:\psi(x)\in\underbrace{[i,i+1)\times[j,j+1)}_{=:\:R_{(i,\:j)}}\right\},$$ i.e. $h=1_B$. Now let $$U_n:=\psi(X_n)\;\;\;\text{for }n\in\mathbb N$$ and note that, by construction of $q$, $$U_n\sim\mathcal U_R\;\;\;\text{for all }n\in\mathbb N\tag1,$$ where $\mathcal U_R$ denotes the uniform distribution on $R$. This immediately yields that $$(q\lambda)(B)=\mathcal U_R\left(R_{(i,\:j)}\right)=\frac1{|J|}\tag2.$$

Let $$g(x):=h(x)\left.\begin{cases}\displaystyle\frac fq(x)&\text{, if }q(x)>0\\0&\text{, otherwise}\end{cases}\right\}\;\;\;\text{for }x\in E.$$ Now the naive estimator for $I$ is $$\frac1n\sum_{i=1}^ng(X_i)\xrightarrow{n\to\infty} I\;\;\;\text{almost surely}\tag3.$$ However, if I implement this (which corresponds to replacing m_camera->film->WriteImage(pbrt::Float{ 1 } / m_sampler->samplesPerPixel); in the last line of importance_sampling_integrator::render() by m_camera->film->WriteImage(pbrt::Float{ 1 } / m_n);), I end up with a black image. This is not surprising, since due to the very small support of $h$ compared to the whole image, only very few terms contribute to the sum in $(3)$.

So, it seems like that one needs to consider only the visits of $(X_n)_{n\in\mathbb N}$ to the set $B$, i.e. consider the process $$Y_k:=X_{\tau_k}\;\;\;\text{for }k\in\mathbb N$$ instead, where $\tau_0:=0$, $$\tau_k:=\inf\left\{n>\tau_{k-1}:X_n\in B\right\}\;\;\;\text{for }k\in\mathbb N.$$ It's easy to see that $(Y_k)_{k\in\mathbb N}$ is independent with $$Y_n\sim (q\lambda)\left[\;\cdot\mid B\right]\;\;\;\text{for all }k\in\mathbb N\tag4.$$ Moreover, $$\frac1k\sum_{i=1}^kg(Y_i)\xrightarrow{k\to\infty}(q\lambda)\left[g\mid B\right]=|J|I\tag5$$ by $(2)$.

Now what the code does is independently drawing $X_1,\ldots,X_n\sim q\lambda$, where $$n=n_0|J|\tag6$$ for some $n_0\in\mathbb N$ (called mutationsPerPixel in PBRT).

We may consider the number $$N_n:=\sum_{i=1}^n1_B(X_i)$$ of visits to $B$ up to time $n$ and note that $$\operatorname E\left[N_n\right]=\frac n{|J|}\tag7$$ and $$\frac{N_n}n\xrightarrow{n\to\infty}\frac1{|J|}\;\;\;\text{almost surely}.\tag8$$

Now PBRT (and the code I've provided above as well) uses the estimator $$\frac1{n_0}\sum_{i=1}^ng(X_i)=\frac1{\operatorname E\left[N_n\right]}\sum_{i=1}^{N_n}g(Y_i),$$ but I absolutely don't understand why this works. From $(5)$ we should obtain $$\frac1{\operatorname E\left[N_n\right]}\sum_{i=1}^{N_n}g(Y_i)=\frac{|J|}n\sum_{i=1}^ng(X_i)\xrightarrow{n\to\infty}|J|I\tag9$$ almost surely, but the left-hand side of $(9)$, dividied by $|J|$, is once again the pracitcally not working estimator $(3)$.

What am I missing?

EDIT: Just to be completely sure I've built the image array by myself using the (supposed to be correct) estimator $(5)$:

void render(pbrt::Scene const& scene)
{
    auto const light_distribution = pbrt::ComputeLightPowerDistribution(scene);

    std::unordered_map<pbrt::Light const*, std::size_t> light_to_index;
    for (std::size_t i = 0; i < scene.lights.size(); ++i)
        light_to_index[scene.lights[i].get()] = i;

    std::vector<vertex> z(m_t),
        y(m_s);

    std::size_t const a = m_camera->film->GetSampleBounds().pMax.y,
        b = m_camera->film->GetSampleBounds().pMax.x;
    std::vector<pbrt::Spectrum> image(a * b);
    std::vector<std::size_t> n(a * b);

    pbrt::MemoryArena arena;
    for (std::size_t i = 0; i < m_n; ++i)
    {
        pbrt::Point2f raster_point;
        auto const contribution = sample_path(scene, *m_camera, *light_distribution, light_to_index,
            arena, m_k, m_s, *m_sampler, raster_point, z.data(), y.data());
        arena.Reset();

        std::size_t const index = static_cast<std::size_t>(raster_point.y) * b + static_cast<std::size_t>(raster_point.x);

        image[index] += contribution.value_or(0);
        ++n[index];
    }

    for (std::size_t i = 0; i < n.size(); ++i)
    {
        if (n[i] > 0)
            image[i] = image[i] / (n[i] /** a * b*/);
    }
    m_camera->film->SetImage(image.data());
    m_camera->film->WriteImage();
}

At the end, I should devide each pixel value by $|J|=ab$. However, if I uncomment the devision in the code, I end up with a black image again. So, my guess right now is that the values provided by BxDF::Sample_f() are not what I thought they are. Maybe they are "normalized" and maybe this is what they mean by "Each Metropolis iteration within <> has splatted contributions with weighted unit luminance".