## 16.4 Metropolis Light Transport

In 1997, Veach and Guibas proposed an unconventional rendering technique named
*Metropolis Light Transport* (MLT), which applies the Metropolis-Hastings
algorithm from Section 13.4 to the path space integral
in Equation (16.1). Whereas all of the rendering techniques we
have discussed
until now have been based on the principles of Monte Carlo integration and
independent sample generation, MLT adopts a different set of tools that allows
samples to be *statistically correlated*.

MLT generates a sequence of light-carrying paths through the scene, where each path is found by mutating the previous path in some manner. These path mutations are done in a way that ensures that the overall distribution of sampled paths is proportional to the contribution the paths make to the image being generated. Such a distribution of paths can in turn be used to generate an image of the scene. Given the flexibility of the Metropolis sampling method, there are relatively few restrictions on the types of admissible mutation rules: highly specialized mutations can be used to sample otherwise difficult-to-explore families of light-carrying paths in a way that would be tricky or impossible to realize without introducing bias in a standard Monte Carlo context.

The correlated nature of MLT provides an important advantage over methods based
on independent sample generation, in that MLT is able to perform a *local
exploration* of path space: when a path that makes a large contribution to the
image is found, it’s easy to find other similar paths by applying small
perturbations to it (a sampling process that generates states based on the
previous state value in this way is referred to as a *Markov chain*). The
resulting short-term memory is often beneficial:
when a function has a small value over most of
its domain and a large contribution in only a small subset, local exploration
amortizes the expense (in samples) of the search for the important region by
taking many samples from this part of the path space. This property makes MLT
a good choice for rendering particularly challenging scenes: while it is
generally not able to match the performance of uncorrelated integrators for
relatively straightforward lighting problems, it distinguishes itself in more
difficult settings where most of the light transport happens along a small
fraction of all of the possible paths through the scene.

The original MLT technique by Veach and Guibas builds on the path space theory of light transport, which presents additional challenges compared to the previous simple examples in Section 13.4.4: path space is generally not a Euclidean domain, as surface vertices are constrained to lie on 2D subsets of . Whenever specular reflection or refraction occurs, a sequence of three adjacent vertices must satisfy a precise geometric relationship, which further reduces the available degrees of freedom.

MLT builds on a set of five mutation rules that each targets specific families of light paths. Three of the mutations perform a localized exploration of particularly challenging path classes such as caustics or paths containing sequences of specular-diffuse-specular interactions, while two perform larger steps with a corresponding lower overall acceptance rate. Implementing the full set of MLT mutations is a significant undertaking: part of the challenge is that none of the mutations is symmetric; hence an additional transition density function must be implemented for each one. Mistakes in any part of the system can cause subtle convergence artifacts that are notoriously difficult to debug.

### 16.4.1 Primary Sample Space MLT

In 2002, Kelemen et al. presented a rendering technique that is also based on
the Metropolis-Hastings algorithm. We will refer to it as *Primary Sample
Space MLT* (PSSMLT) for reasons that will become clear shortly. Like MLT, the
PSSMLT method explores the space of light paths, searching for
those that carry a significant amount of energy from a light source to the
camera. The main difference is that PSSMLT does not use path space directly: it
explores light paths indirectly, piggybacking on an existing Monte Carlo
rendering algorithm such as uni- or bidirectional path tracing, which has both
advantages and disadvantages. The main advantage is that PSSMLT can use
symmetric transitions on a Euclidean domain, and for this reason it is much
simpler to implement. The downside is that PSSMLT lacks detailed information
about the structure of the constructed light paths, making it impossible to
recreate the kinds of sophisticated mutation strategies that are found in the
original MLT method.

The details of this approach are easiest to motivate from an
implementation-centric viewpoint. Consider the case of rendering a scene
with the path integrator from Section 14.5.4, where
we generate (pseudo-)random samples in raster space to get a starting point
on the film, turn
them into rays, and invoke `PathIntegrator::Li()` to obtain
corresponding radiance estimates.
`PathIntegrator::Li()` will generally request a number of additional 1D
or 2D samples from the `Sampler` to sample light sources and
material models, perform Russian roulette tests, and so on. Conceptually,
these samples could all be generated ahead of time and passed to
`PathIntegrator::Li()` using an additional argument, thus completely
removing any (pseudo-)randomness from an otherwise purely deterministic
function.

Suppose that is the resulting deterministic function, which maps an infinite-dimensional sample sequence to a radiance estimate . Here, denotes the raster position, and the remaining arguments are the samples consumed by . By drawing many samples and averaging the results of multiple evaluations of , we essentially compute high-dimensional integrals of over a “hypercube” of (pseudo-)random numbers

Given this interpretation of integrating a radiance estimate function over
a domain made of all possible sequences of input samples, we can apply the
Metropolis-Hastings algorithm to create a sampling process on the same
domain, with a distribution that is proportional to
(Figure 16.19). This distribution is intuitively
desirable: more sampling effort is naturally placed in parts of the
sampling space where more light transport occurs. The state space
for this problem consists of infinite dimensional sample
vectors and is called
the *primary sample space*. For brevity in the equations below, we
will write the state vector as

We use the convention that all components (including the raster coordinates) are in the interval and appropriately re-scaled within when necessary.

PSSMLT explores primary sample space using two different kinds of mutations. The first (the “large step” mutation) replaces all components of the vector with new uniformly distributed samples, which corresponds to invoking the underlying Monte Carlo rendering method as usual (i.e., without PSSMLT).

Recall from Section 13.4.2 that it is important that there be greater than zero probability that any possible sample value be proposed; this is taken care of by the large step mutation. In general, the large step mutation helps us explore the entire state space without getting stuck on local “islands.”

The second mutation (the “small step” mutation) makes a small perturbation to each of the sample values . This mutation explores light-carrying paths close to the current path, which is particularly important when a difficult lighting configuration is encountered.

Both of these are symmetric mutations, so their transition probabilities cancel out when the acceptance probability is computed and thus don’t need to be computed, as was shown in Equation (13.8).

The interface between the outer Metropolis-Hastings iteration and the inner
`Integrator` only involves the exchange of abstract sample vectors, making
this an extremely general approach: PSSMLT can theoretically enhance any kind
of rendering method that is based on Monte Carlo integration—in fact, it even
works for general Monte Carlo integration problems that are not related to
rendering at all.

In practice, PSSMLT is often implemented on top of an existing bidirectional path tracer. The resulting algorithm generates a new primary sample space state in every iteration and passes it to BDPT, which invokes its set of connection strategies and re-weights the result using MIS. In this setting, mutations effectively jump from one group of path connections to another group rather than dealing with individual connection strategies. However, this is not without disadvantages: in many cases, only a small subset of the strategies is truly effective, and MIS will consequently assign a large weight only to this subset. This means that the algorithm still spends a considerable portion of its time generating connections with strategies that have a low weight and thus contribute very little to the rendered image.

### 16.4.2 Multiplexed MLT

In 2014, Hachisuka et al. presented an extension to PSSMLT named
*Multiplexed Metropolis Light Transport* (MMLT) to address this problem.
MMLT leaves the “outer” Metropolis-Hastings iteration conceptually unchanged
and applies a small but effective modification to the “inner” BDPT
integrator. Instead of always invoking all BDPT connection strategies, the
algorithm chooses a single strategy according to an extra state dimension and
returns its contribution scaled by the inverse discrete probability of the
choice. The additional dimension used for strategy selection can be mutated
using small or large steps in the same way as the other primary sample space
components of .

To prevent unintentional large structural path mutations, Hachisuka et al. fix the Markov chain so that it only explores paths of a fixed depth value; the general light transport problem is then handled by running many independent Markov chains.

The practical consequence is that the Metropolis sampler will tend to spend more computation on effective strategies that produce larger MIS-weighted contributions to the image. Furthermore, the individual iterations are much faster since they only involve a single connection strategy. The combination of these two aspects improves the Monte Carlo efficiency of the resulting estimator.

Figure 16.20 shows the contemporary house scene rendered with bidirectional path tracing and MMLT, using roughly equal computation time for each approach.

*(Model courtesy of Florent Boyer.)*

Figure 16.21 compares path tracing, BDPT, and MMLT with the San Miguel scene. For both scenes, MMLT generates a better result, but the difference is more pronounced with the house scene, which is particularly challenging for light transport algorithms in that there is essentially no direct illumination inside the house; all light-carrying paths must follow specular bounces through the glass windows.

*(Model courtesy of Guillermo M. Leal Llaguno.)*

Table 16.1 helps explain these efficiency differences: for both of these scenes, both path tracing and BDPT have a lot of trouble finding paths that carry any radiance, while Metropolis is more effective at doing so thanks to path reuse.

Path tracing | BDPT | MMLT | |
---|---|---|---|

Modern House | 98.0% | 97.6% | 51.9% |

San Miguel | 95.9% | 97.0% | 62.0% |

### 16.4.3 Application to Rendering

Metropolis sampling generates samples from the distribution of a given scalar function. To apply it to rendering, there are two issues that we must address: first, we need to estimate a separate integral for each pixel to turn the generated samples into an image, and, second, we need to handle the fact that is a spectrally valued function but Metropolis needs a scalar function to compute the acceptance probability, Equation (13.7).

We can apply the ideas from Section 13.4.5 to use Metropolis
samples to compute integrals. First, we define the *image contribution
function* such that for an image with pixels, each pixel has a value
that is the integral of the product of the pixel’s image reconstruction filter
and the radiance that contributes to the image:

The filter function only depends on the two components of associated with the raster position. The value of for any particular pixel is usually 0 for the vast majority of samples due to the filter’s finite extent.

If samples are generated from some distribution, , then the standard Monte Carlo estimate of is

Recall that Metropolis sampling requires a scalar function that defines the
desired distribution of samples generated by the algorithm. Unfortunately,
is a spectrally valued function, and thus there is no unambiguous notion of
what it means to generate samples proportional to . To work around this
issue, we will define a *scalar contribution function* that
is used inside the Metropolis iteration. It is desirable that this function be
large when is large so that the distribution of samples has some
relationship to the important regions of . As such, using the luminance of
the radiance value is a good choice for the scalar contribution function. In
general, any function that is nonzero when is nonzero will generate correct
results, just possibly not as efficiently as a function that is more directly
proportional to .

Given a suitable scalar contribution function, , Metropolis generates a sequence of samples from ’s distribution, the normalized version of :

and the pixel values can thus be computed as

The integral of over the entire domain can be computed using a traditional approach like bidirectional path tracing. If this value is denoted by , with , then each pixel’s value is given by

In other words, we can use Metropolis sampling to generate samples from the distribution of the scalar contribution function . For each sample, the pixels it contributes to (based on the extent of the pixel filter function ) have the value

added to them. Thus, brighter pixels have larger values than dimmer pixels due to more samples contributing to them (as long as the ratio is generally of the same magnitude).

### 16.4.4 Primary Sample Space Sampler

The `MLTIntegrator` applies Metropolis sampling and MMLT to render
images, using the bidirectional path tracer from
Section 16.3. Its implementation is contained in
the files `integrators/mlt.h` and
`integrators/mlt.cpp`.
Figure 16.22 shows the effectiveness of this integrator
with a particularly tricky lighting situation.

`MLTIntegrator`, using a roughly equal amount of time. Note that MMLT gives a lower variance result, thanks to being able to efficiently explore the local path space once a high-contribution path has been found.

Before describing the `MLTIntegrator` implementation, we’ll first introduce the `MLTSampler`,
which is responsible for managing primary sample space state vectors,
mutations, and acceptance and rejection steps.

The `MLTIntegrator` works best if the `MLTSampler` actually
maintains three separate sample vectors—one for the camera subpath, one
for the light subpath, and one for the connection step. We’ll say that
these are three *sample streams*. The
`streamCount` parameter to the constructor
lets the caller request a particular number of
such sample streams.

Later, during the initialization phase of `MLTIntegrator`, we will
create many separate `MLTSampler` instances that are used to select a
suitable set of starting states for the Metropolis sampler. Importantly,
this process requires that each `MLTSampler` produces a distinct
sequence of state vectors. The `RNG` pseudo-random number generator
used by `pbrt` has a handy feature that makes this easy to accomplish: the
`RNG` constructor accepts a *sequence index*, which selects between
one of unique pseudo-random sequences. We thus add an
`rngSequenceIndex` parameter to the `MLTSampler` constructor that
is used to supply a unique stream index to the internal `RNG`.

The `largeStepProbability` parameter refers to the probability of taking a
“large step” mutation, and `sigma` controls the size of “small step”
mutations.

The `MLTSampler::X` member variable stores the current sample vector
. Because we generally don’t know ahead of time how many
dimensions of are actually needed during the sampler’s lifetime,
we’ll start with an empty vector and expand it on demand as calls to
`MLTSampler::Get1D()` and `MLTSampler::Get2D()` occur during
rendering.

The elements of this array have the type `PrimarySample`. The main task for
`PrimarySample` is to record the current value of a single component of
on the interval . In the following, we’ll add some
additional functionality for representing proposed mutations and restoring
the original sample value if a proposed mutation is rejected.

The `Get1D()` method returns the value of a single component of
`MLTSampler::X`, whose position is given by
`GetNextIndex()`—for now, we can think of this method as
returning the value of a running counter that increases with every call.
`EnsureReady()` expands `MLTSampler::X` as needed and
ensures that its contents are in a consistent state. The details of these
methods will be clearer after a few more preliminaries, so we won’t
introduce their implementations just yet.

The 2D analog simply performs two calls to `Get1D()`.

Next, we will define several `MLTSampler` methods that
are not part of the official `Sampler` interface. The first one,
`MLTSampler::StartIteration()`, is called at the beginning of each
Metropolis-Hastings iteration; it increases the `currentIteration`
counter and determines which type of mutation (small or large) should be
applied to the sample vector in the current iteration.

`currentIteration` is a running counter, which keeps track of the
current Metropolis-Hastings iteration index. Note that iterations with
rejected proposals will be excluded from this count.

The `MLTSampler::lastLargeStepIteration` member variable records the
index of the last iteration where a successful large step took place. Our
implementation chooses the initial state to be uniformly
distributed on primary sample space ; hence the first iteration’s
state can be interpreted as the result of a large step in iteration 0.

The `MLTSampler::Accept()` method is called whenever a proposal is
accepted.

At this point, the main missing parts are `MLTSampler::EnsureReady()` and
the logic that applies the actual mutations to `MLTSampler::X`.
Before filling in those gaps, let us take a brief step back.

In theory, all entries in the `X` vector must be updated by small or
large mutations in every iteration of the Metropolis sampler. However,
doing so would sometimes be rather inefficient: consider the case where most
iterations only query a small number of dimensions of , hence
`MLTSampler::X` has not grown to a large size yet. If a later iteration
makes many calls to `Get1D()` or `Get2D()`,
then the dynamic array `MLTSampler::X` must correspondingly expand, and
these extra entries increase the cost of every subsequent Metropolis
iteration (even if these components of are never referenced
again!).

Instead, it’s more efficient to update the `PrimarySample` entries on
demand—that is, when they are referenced by an actual `Get1D()` or
`Get2D()` call. Doing so avoids the aforementioned inefficiency,
though after some period of inactivity, we must carefully replay all
mutations that a given `PrimarySample` missed. To keep track of this
information, an additional member variable recording the last iteration
where a `PrimarySample` was modified is useful.

We now have enough background to proceed to the implementation of the
`MLTSampler::EnsureReady()` method, which updates individual sample values
when they are accessed.

`MLTSampler::X`if necessary and get current >> <<Reset if a large step took place in the meantime>>

`sample`>>

`nSmall`small step mutations>>

First, any gaps are filled with zero-initialized `PrimarySample`s before
a reference to the requested entry is obtained.

`MLTSampler::X`if necessary and get current >>=

When the last modification of `Xi` precedes the last large step,
the current content of `Xi.value` is irrelevant, since it should
have been overwritten with a new uniform sample in iteration
`lastLargeStepIteration`. In this case, we simply replay this missed
mutation and update the last modification iteration index accordingly.

Next, a call to `Backup()` notifies the `PrimarySample` that a
mutation is going to be proposed; doing so allows it to make a copy of
’s sample value in case the mutation is rejected. All remaining
mutations between iterations `lastLargeStepIteration` and
`currentIteration` are then applied. Two different cases can occur:
when the current iteration is a large step, we simply initialize
`PrimarySample::value` with a uniform sample. Otherwise, all iterations
since the last large step are (by definition) small steps that we must
replay.

`sample`>>=

`nSmall`small step mutations>>

For small steps, we apply normally distributed perturbations to each component:

where is given by `MLTIntegrator::sigma`.
The advantage of sampling with a normal distribution like this is that it
naturally tries a variety of mutation sizes. It preferentially makes small
mutations that remain close to the current state, which help locally explore
the path space in small areas of high contribution where large mutations would
tend to be rejected. On the other hand, because it also can make larger
mutations, it also avoids spending too much time in a small part of the path
space in cases where larger mutations have a good likelihood of acceptance.

Our implementation here merges the sequence of `nSmall`
perturbations into a single update and clamps the result to the unit interval,
wrapping around to the other end of the domain if necessary. Wrapping ensures
that the transition probabilities for all pairs of sample values are symmetric.

`nSmall`small step mutations>>=

Recall the rule that when two samples from normal distributions and are added, the sum is also normally distributed with parameters . Thus, when perturbations are to be applied to , instead of performing perturbations in sequence, it is equivalent and more efficient to directly sample

which we do by importance sampling a standard normal distribution and scaling the result by . Applying the inversion method to the PDF

gives the following sampling method for a uniform sample :

where is the error function, , and is its inverse.
The `ErfInv()` function, not included here, approximates
with a polynomial.

We scale the resulting sample by the effective variance from Equation (16.23) and use it to perturb the sample before keeping only its fractional component, so that it remains in .

`MLTSampler::Reject()` must be called whenever a proposed mutation is
rejected. It restores all `PrimarySample`s modified in the current
iteration and reverts the iteration counter.

The `Backup()` and `Restore()` methods make it possible to record
the value of a `PrimarySample` before a mutation and to restore it if
the mutation is rejected.

Before wrapping up `MLTSampler`, we must address a detail that would
otherwise cause issues when the sampler is used with BDPT. For each pixel
sample, the `BDPTIntegrator` implementation calls
`GenerateCameraSubpath()` and `GenerateLightSubpath()` in sequence
to generate a pair of subpaths, each function requesting 1D and 2D samples
from a supplied `Sampler` as needed.

In the context of MLT, the resulting sequence of sample requests creates a mapping between components of and vertices on the camera or light subpath. With the process described above, the components determine the camera subpath (for some ), and the remaining values determine the light subpath. If the camera subpath requires a different number of samples after a perturbation (e.g., because the random walk produced fewer vertices), then there is a shift in the assignment of primary sample space components to the light subpath. This leads to an unintended large-scale modification to the light path.

It is easy to avoid this problem with a more careful indexing scheme: the
`MLTSampler` partitions into multiple interleaved streams
that cannot interfere with each other.
The `MLTSampler::StartStream()` method indicates that subsequent
samples should come from the stream with the given index. It also resets
`sampleIndex`, the index of the current sample in the stream. (The
number of such streams, `MLTSampler::streamCount`, is specified in the
`MLTSampler` constructor.)

After the stream is selected, the `MLTSampler::GetNextIndex()` method
performs corresponding steps through the primary sample vector components.
It interleaves the streams into the global sample vector—in other words,
the first `streamCount` dimensions in are respectively used
for the first dimension of each of the streams, and so forth.

### 16.4.5 MLT Integrator

Given all of this infrastructure—an explicit representation of an
-dimensional sample , functions to apply mutations to it, and
BDPT’s vertex abstraction layer for evaluating the radiance of a given sample
value—we can move forward to the heart of the implementation of the
`MLTIntegrator`.

The `MLTIntegrator` constructor, not shown here, just initializes
various member variables from parameters provided to it. These member
variables will be introduced in the following as they are used.

We begin by defining the method `MLTIntegrator::L()`, which computes
the radiance for a vector of sample values provided
by an `MLTSampler`. Its parameter `depth` specifies a specific path
depth, and `pRaster` returns the raster position of the path, if the
path successfully carries light from a light source to the film plane. The
initial statement activates the first of three streams in the underlying
`MLTSampler`.

`t`vertices>>

`s`vertices>>

The implementation uses three sample streams from the `MLTSampler`: the
first two for the camera and light subpath and the third one for any
`Camera::Sample_Wi()` or `Light::Sample_Li()` calls performed by
connection strategies in `ConnectBDPT()` with or (refer to
Section 16.3.3 for details).

The body of `MLTIntegrator::L()` first selects an individual BDPT
strategy for the provided `depth` value—this is the MMLT
modification to PSSMLT—and invokes the bidirectional path-tracing
machinery to compute a corresponding radiance estimate. For paths with zero
scattering events (i.e., directly observed light sources), the only viable
strategy provided by the underlying BDPT implementation entails
intersecting them with a ray traced from the camera. For
longer paths, there are possible strategies. The fragment
below uniformly maps the first primary sample space dimension onto this
set of strategies. The variables `s` and `t` denote the number of
light and camera subpath sampling steps following the convention of the
`BDPTIntegrator`.

The next three fragments compute the radiance estimate. They strongly
resemble BDPT with some MMLT-specific modifications: the first one samples
a film position in , maps it to raster coordinates, and tries to
generate a corresponding camera subpath with *exactly* vertices,
failing with a 0-valued estimate for when this was not possible.

`t`vertices>>=

The `camera` member variable holds the `Camera` specified in the
scene description file.

The next fragment implements an analogous operation for the light
subpath. Note the call to `MLTSampler::StartStream()`, which switches
to the second stream of samples.

`s`vertices>>=

Finally, we switch to the last sample stream and invoke the strategy
via a call to `ConnectBDPT()`. The final radiance estimate is multiplied by the
inverse probability of choosing the current strategy , which is
equal to `nStrategies`.

#### Main Rendering Loop

There are two phases to the rendering process implemented in
`MLTIntegrator::Render()`. The first phase generates a set of bootstrap
samples that are candidates for initial states of Markov chains and
computes the normalization constant ,
from Equation (16.22). The second phase runs a series
of Markov chains, where each chain chooses one of the bootstrap samples for
its initial sample vector and then applies Metropolis sampling.

`i`th bootstrap sample>>

`nChains`Markov chains in parallel>>

`nChainMutations`>>

`nChainMutations`steps>>

`film`>>

Following the approach described in Section 13.4.3, we avoid issues with start-up bias by computing a set of bootstrap samples with a standard Monte Carlo estimator and using them to create a distribution that supplies the initial states of the Markov chains. This process builds on the sample generation and evaluation routines implemented previously.

Each bootstrap sample is technically a sequence of
samples with different path depths; the following fragment initializes the
array `bootstrapWeights` with their corresponding luminance values. At
the end, we create a `Distribution1D` instance to sample bootstrap
paths proportional to their luminances and set the constant `b` to the
sum of average luminances for each depth. Because we’ve kept the
contributions of different path sample depths distinct, we can
preferentially sample path lengths that make the largest contribution to
the image. It is straightforward to compute the bootstrap samples in
parallel, as all loop iterations are independent from one another.

`i`th bootstrap sample>>

As usual, `maxDepth` denotes the maximum number of interreflections that
should be considered. `nBootstrap` sets the number of
bootstrapping samples to use to seed the iterations and compute the integral
of the scalar contribution function, Equation (16.22).

In each iteration, we instantiate a dedicated `MLTSampler` with index
`rngIndex` providing a unique sample vector that is
uniformly distributed in the primary sample space. Next, we evaluate to
obtain a corresponding radiance estimate for the current path depth
`depth` and write its luminance into `bootstrapWeights`. The
raster-space position `pRaster` is not needed in the preprocessing
phase, so it is ignored.

`i`th bootstrap sample>>=

The `mutationsPerPixel` parameter is analogous to
`Sampler::samplesPerPixel` and denotes the number of iterations that MLT
(on average!) spends in each pixel. Individual pixels will receive an actual
number of samples related to their brightness, due to Metropolis’s property of
taking more samples in regions where the function’s value is high. The total
number of Metropolis samples taken is the product of the number of pixels
and `mutationsPerPixel`.

The `sigma` and `largeStepProbability` member variables give the
corresponding configuration parameters of the `MLTSampler`.

We now move on to the main rendering task, which performs a total of
`nTotalMutations` mutation steps spread out over `nChains` parallel
Markov chains.

We must be careful that the actual number of mutation steps indeed
comes out to be equal to `nTotalMutations`, particularly when `nTotalMutations`
is not divisible by the number of parallel chains. The solution is simple: we
potentially stop the last Markov chain a few iterations short; the corrected
number of per-chain iterations is given by `nChainMutations`.

`nChains`Markov chains in parallel>>=

`nChainMutations`>>

`nChainMutations`steps>>

`film`>>

`nChains` specifies the number of Markov chains that should be executed
independently of each other. Its default value of 100 is a trade-off between
providing sufficient parallelism and running each chain for a long time.

The MLT integrator only splats contributions to arbitrary pixels in the
film; it doesn’t fill in well-defined tiles of the image plane. Therefore,
a `FilmTile` isn’t necessary, and calls to `Film::AddSplat()`
suffice to update the image.

`nChainMutations`>>=

`nChainMutations`steps>>

`film`>>

Every Markov chain instantiates its own pseudo-random number generator
following a unique stream. Note that this `RNG` instance is separate from
the one in `MLTSampler`: it is used to pick the initial state and to
accept or reject Metropolis proposals later on. Due to the ordering used
earlier when initializing the
entries of `bootstrap` array, we can immediately deduce the path depth of the
sampled index `bootstrapIndex`.
An important consequence of the method used to generate the `bootstrap`
distribution is that the expected number of sampled initial states with a given
depth value is proportional to their contribution to the image.

Having chosen the bootstrap sample, we must now obtain the corresponding
primary sample space vector . One way to accomplish this would
have been to store all `MLTSampler` instances in the preprocessing
phase. A more efficient approach builds on the property that the
initialization in the `MLTSampler` constructor is completely
deterministic and just depends on the `rngSequenceIndex`
parameter. Thus, here we can create an `MLTSampler` with index
`bootstrapIndex`, which recreates the same sampler that
originally produced the sampled entry of the `bootstrap` distribution.

With the sampler in hand, we can compute the current value of and its position on the film.

The implementation of the Metropolis sampling routine in the following fragments follows the expected values technique from Section 13.4.1: a mutation is proposed, the value of the function for the mutated sample and the acceptance probability are computed, and the weighted contributions of both the new and old samples are recorded. The proposed mutation is then randomly accepted based on its acceptance probability.

One of the unique characteristics of MMLT
(Section 16.4.2) compared to
other Metropolis-type methods is that each Markov chain is restricted to paths
of a fixed `depth` value. The first sample dimension selects among
various different strategies, but only those producing equal path depths are
considered, which improves performance of the method by making proposals more
local. The contribution of all path depths is accounted for by starting many
Markov chains with different initial states.

`nChainMutations`steps>>=

`film`>>

Given the scalar contribution function’s value, the acceptance probability is then given by the simplified expression from Equation (13.8) thanks to the symmetry of our mutations on primary sample space.

As described at the start of the section, the spectral radiance value must be converted to a value given by the scalar contribution function so that the acceptance probability can be computed for the Metropolis sampling algorithm. Here, we compute the path’s luminance, which is a reasonable choice.

Both samples can now be added to the image. Here, they are scaled with weights based on the expected values optimization introduced in Section 13.4.1.

`film`>>=

Finally, the proposed mutation is either accepted or rejected, based on the
computed acceptance probability `accept`. If the mutation is accepted,
then the values `pProposed` and `LProposed` become properties of
the current state. In either case, the `MLTSampler` must be informed of the
outcome so that it can update the `PrimarySample` accordingly.

Metropolis sampling only considers the *relative* frequency of samples
and cannot create an image that is correctly scaled in *absolute*
terms; hence the value is crucial: it contains an estimate of the
average luminance of the `Film` that we use to remove this
ambiguity. 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`. We
thus must cancel this factor out and multiply by when writing the image
to convert to actual incident radiance on the film.