Return to Question

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?

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?

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();
}

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();
}

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".