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.

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.

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.

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.

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.

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

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.

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.

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

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.

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.

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.

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.

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.

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

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.

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. ^†

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.