# Scaling of the final image in Metropolis Light Transport

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

• Are you certain they have what you refer to as $g$ in their code/formulation? It could be that they have $\frac{hf}{q'}$ where $q'$ is the probability without the first vertex (for which it is $\frac{1}{|J|}$). If that is the case then what you wrote and what they have would be equivalent, since what you have is really: $\frac{1}{spp}\sum \frac{h(X_i)f(X_i)}{q'(X_i)}$. $q'$ would be the probability of picking the remaining dimensions of the sample (without the factor coming from picking a uniformly a sample on the film) - that is $q(X_i) = \frac{q'(X_i)}{|J|}$. Mar 5 '20 at 10:26
• With regards to your new edit - how does sample_path compute the probability? If it doesn't include the $\frac{1}{|J|}$ probability multiplication (the contribution being divided by it), then what happens is exactly what I outlined above. Mar 5 '20 at 10:47
• @lightxbulb That's a good point. Well, they don't compute $f$ and $q$ separately, but only the fraction $\frac fq$ (which is what I called contribution and they refer to a $\beta$ in the book). The relevant $β_1⋯β_k$ here would be the one for the second vertex $z_1$ on the camera subpath. This $\beta$ is always equal to $1$. Formally, it should be $$\beta_1:=\frac{W_e(z_1\to z_0)}{p_0(z_0,\omega_{z_0\to z_1})},$$ where $p_0(z_0,\omega_{z_0\to z_1})$ is the density of choosing the first ray $r_1=(z_0,\omega_{z_0\to z_1})$ on the camera subpath with respect to the path throughput measure. Mar 5 '20 at 10:48
• @lightxbulb Now, if I understood them correctly, they defined $W_e(z_0\to z_1)$ (for a perspective camera) to be equal to this density; see equation 16.3 here: pbr-book.org/3ed-2018/Light_Transport_III_Bidirectional_Methods/…). Maybe the normalization of $W_e(z_0\to z_1)$ is the source of the error. What do you think? Mar 5 '20 at 10:49
• Then they simply do not have this $\frac{1}{|J|}$ in the probability, since they have taken it outside. This means that their contribution is divided by $q'$, instead of $q$, which explains the $\frac{1}{spp}$. They basically simply took out the raster probability out of the sum, and canceled it out. Note that this makes sense, especially since this makes it consistent with the other pixel estimators. In general if you had pixels with varying area, you would instead have $\frac{A_j}{\sum_k A_k}$ for example. Simply put - you're missing the $\frac{1}{|J|}$ in your probability in your code. Mar 5 '20 at 10:53