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 bold upper R cubed . 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 upper L is the resulting deterministic function, which maps an infinite-dimensional sample sequence bold upper X 1 comma bold upper X 2 comma ellipsis to a radiance estimate upper L left-parenthesis bold upper X 1 comma bold upper X 2 comma ellipsis right-parenthesis . Here, left-parenthesis bold upper X 1 comma bold upper X 2 right-parenthesis denotes the raster position, and the remaining arguments bold upper X 3 comma bold upper X 4 comma ellipsis are the samples consumed by upper L . By drawing many samples and averaging the results of multiple evaluations of upper L , we essentially compute high-dimensional integrals of upper L over a “hypercube” of (pseudo-)random numbers bold upper X 1 comma bold upper X 2 comma ellipsis period

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 upper L (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 normal upper Omega for this problem consists of infinite dimensional sample vectors left-parenthesis bold upper X 1 comma bold upper X 2 comma ellipsis right-parenthesis element-of left-bracket 0 comma 1 right-parenthesis Superscript normal infinity and is called the primary sample space. For brevity in the equations below, we will write the state vector as

bold upper X equals left-parenthesis bold upper X 1 comma bold upper X 2 comma ellipsis right-parenthesis period

We use the convention that all components (including the raster coordinates) are in the interval left-bracket 0 comma 1 right-parenthesis and appropriately re-scaled within upper L when necessary.

Figure 16.19: Primary Sample Space MLT. PSSMLT performs mutations in an abstract space of infinite dimensional “random number vectors” bold upper X . A deterministic mapping upper L constructs corresponding light paths on path space and estimates their radiance.

PSSMLT explores primary sample space using two different kinds of mutations. The first (the “large step” mutation) replaces all components of the vector bold upper X 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 bold upper X Subscript i . 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 bold upper X .

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.20: Comparison of Bidirectional Path Tracing and Multiplexed Metropolis Light Transport. (1) Rendered with bidirectional path tracing with 128 samples per pixel. Even with many samples, the image is quite noisy. (2) Rendered with Multiplexed Metropolis Light Transport with an average of 420 mutations per pixel (roughly equal running time). MLT generates a substantially better image for the same amount of work. BDPT does do better in a few localized parts of the image: note the directly illuminated surfaces outside the windows, for example. There, the ability to use well-distributed sample points for the direct lighting calculation gives a result with lower variance. (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.

Figure 16.21: Comparison of Path Tracing, BDPT, and Multiplexed Metropolis Light Transport. (1) Path tracing with 200 samples per pixel. (2) Bidirectional path tracing with 128 samples per pixel. (3) Metropolis Light Transport with an average of 950 mutations per pixel (roughly equal running time). BDPT is notably more effective than path tracing, and as with Figure 16.20, MLT is the most effective. (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.

Table 16.1: Percentage of Traced Paths That Carried Zero Radiance. With these two scenes, both path tracing and BDPT have trouble finding light-carrying paths: the vast majority of the generated paths don’t carry any radiance at all. Thanks to local exploration, Metropolis is better able to find additional light-carrying paths after one has been found. This is one of the reasons why it is more efficient than those approaches here.

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 upper L 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 j pixels, each pixel upper I Subscript j has a value that is the integral of the product of the pixel’s image reconstruction filter h Subscript j and the radiance upper L that contributes to the image:

upper I Subscript j Baseline equals integral Underscript normal upper Omega Endscripts h Subscript j Baseline left-parenthesis bold upper X right-parenthesis upper L left-parenthesis bold upper X right-parenthesis normal d normal upper Omega period

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

If upper N samples bold upper X Subscript i are generated from some distribution, bold upper X Subscript i Baseline tilde p left-parenthesis bold upper X right-parenthesis , then the standard Monte Carlo estimate of upper I Subscript j is

upper I Subscript j Baseline almost-equals StartFraction 1 Over upper N EndFraction sigma-summation Underscript i equals 1 Overscript upper N Endscripts StartFraction h Subscript j Baseline left-parenthesis bold upper X Subscript i Baseline right-parenthesis upper L left-parenthesis bold upper X Subscript i Baseline right-parenthesis Over p left-parenthesis bold upper X Subscript i Baseline right-parenthesis EndFraction period

Recall that Metropolis sampling requires a scalar function that defines the desired distribution of samples generated by the algorithm. Unfortunately, upper L is a spectrally valued function, and thus there is no unambiguous notion of what it means to generate samples proportional to upper L . To work around this issue, we will define a scalar contribution function upper I left-parenthesis bold upper X right-parenthesis that is used inside the Metropolis iteration. It is desirable that this function be large when upper L is large so that the distribution of samples has some relationship to the important regions of upper L . 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 upper L is nonzero will generate correct results, just possibly not as efficiently as a function that is more directly proportional to upper L .

Given a suitable scalar contribution function, upper I left-parenthesis bold upper X right-parenthesis , Metropolis generates a sequence of samples bold upper X Subscript i from upper I ’s distribution, the normalized version of upper I :

p left-parenthesis bold upper X right-parenthesis equals StartFraction upper I left-parenthesis bold upper X right-parenthesis Over integral Underscript normal upper Omega Endscripts upper I left-parenthesis bold upper X right-parenthesis normal d normal upper Omega EndFraction comma

and the pixel values can thus be computed as

upper I Subscript j Baseline almost-equals StartFraction 1 Over upper N EndFraction sigma-summation Underscript i equals 1 Overscript upper N Endscripts StartFraction h Subscript j Baseline left-parenthesis bold upper X Subscript i Baseline right-parenthesis upper L left-parenthesis bold upper X Subscript i Baseline right-parenthesis Over upper I left-parenthesis bold upper X Subscript i Baseline right-parenthesis EndFraction left-parenthesis integral Underscript normal upper Omega Endscripts upper I left-parenthesis bold upper X right-parenthesis normal d normal upper Omega right-parenthesis period

The integral of upper I over the entire domain normal upper Omega can be computed using a traditional approach like bidirectional path tracing. If this value is denoted by b , with b equals integral upper I left-parenthesis bold upper X right-parenthesis normal d normal upper Omega , then each pixel’s value is given by

upper I Subscript j Baseline almost-equals StartFraction b Over upper N EndFraction sigma-summation Underscript i equals 1 Overscript upper N Endscripts StartFraction h Subscript j Baseline left-parenthesis bold upper X Subscript i Baseline right-parenthesis upper L left-parenthesis bold upper X Subscript i Baseline right-parenthesis Over upper I left-parenthesis bold upper X Subscript i Baseline right-parenthesis EndFraction period

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

StartFraction b Over upper N EndFraction StartFraction h Subscript j Baseline left-parenthesis bold upper X Subscript i Baseline right-parenthesis upper L left-parenthesis bold upper X Subscript i Baseline right-parenthesis Over upper I left-parenthesis bold upper X Subscript i Baseline right-parenthesis EndFraction

added to them. Thus, brighter pixels have larger values than dimmer pixels due to more samples contributing to them (as long as the ratio upper L Subscript i Baseline left-parenthesis bold upper X Subscript i Baseline right-parenthesis slash upper I left-parenthesis bold upper X Subscript i Baseline right-parenthesis 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.

Figure 16.22: A Volumetric Caustic. Light passing through the sphere is focused in the medium behind it, creating a volumetric caustic. (a) Rendered with bidirectional path tracing, (b) rendered with the 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.

<<MLTSampler Declarations>>= 
class MLTSampler : public Sampler { public: <<MLTSampler Public Methods>> 
MLTSampler(int mutationsPerPixel, int rngSequenceIndex, Float sigma, Float largeStepProbability, int streamCount) : Sampler(mutationsPerPixel), rng(rngSequenceIndex), sigma(sigma), largeStepProbability(largeStepProbability), streamCount(streamCount) { } Float Get1D(); Point2f Get2D(); std::unique_ptr<Sampler> Clone(int seed); void StartIteration(); void Accept(); void Reject(); void StartStream(int index); int GetNextIndex() { return streamIndex + streamCount * sampleIndex++; }
protected: <<MLTSampler Private Declarations>> 
struct PrimarySample { Float value = 0; <<PrimarySample Public Methods>>  <<PrimarySample Public Data>> 
int64_t lastModificationIteration = 0; Float valueBackup = 0; int64_t modifyBackup = 0;
};
<<MLTSampler Private Methods>> 
void EnsureReady(int index);
<<MLTSampler Private Data>> 
RNG rng; const Float sigma, largeStepProbability; const int streamCount; std::vector<PrimarySample> X; int64_t currentIteration = 0; bool largeStep = true; int64_t lastLargeStepIteration = 0; int streamIndex, sampleIndex;
};

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 2 Superscript 63 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.

<<MLTSampler Public Methods>>= 
MLTSampler(int mutationsPerPixel, int rngSequenceIndex, Float sigma, Float largeStepProbability, int streamCount) : Sampler(mutationsPerPixel), rng(rngSequenceIndex), sigma(sigma), largeStepProbability(largeStepProbability), streamCount(streamCount) { }

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

<<MLTSampler Private Data>>= 
RNG rng; const Float sigma, largeStepProbability; const int streamCount;

The MLTSampler::X member variable stores the current sample vector bold upper X . Because we generally don’t know ahead of time how many dimensions of bold upper X 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.

<<MLTSampler Private Data>>+=  
std::vector<PrimarySample> X;

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 bold upper X on the interval left-bracket 0 comma 1 right-parenthesis . 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.

<<MLTSampler Private Declarations>>= 
struct PrimarySample { Float value = 0; <<PrimarySample Public Methods>>  <<PrimarySample Public Data>> 
int64_t lastModificationIteration = 0; Float valueBackup = 0; int64_t modifyBackup = 0;
};

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.

<<MLTSampler Method Definitions>>= 
Float MLTSampler::Get1D() { int index = GetNextIndex(); EnsureReady(index); return X[index].value; }

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

<<MLTSampler Method Definitions>>+=  
Point2f MLTSampler::Get2D() { return Point2f(Get1D(), 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.

<<MLTSampler Method Definitions>>+=  
void MLTSampler::StartIteration() { currentIteration++; largeStep = rng.UniformFloat() < largeStepProbability; }

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.

<<MLTSampler Private Data>>+=  
int64_t currentIteration = 0; bool largeStep = true;

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 bold upper X 0 to be uniformly distributed on primary sample space left-bracket 0 comma 1 right-parenthesis Superscript normal infinity ; hence the first iteration’s state can be interpreted as the result of a large step in iteration 0.

<<MLTSampler Private Data>>+=  
int64_t lastLargeStepIteration = 0;

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

<<MLTSampler Method Definitions>>+=  
void MLTSampler::Accept() { if (largeStep) lastLargeStepIteration = currentIteration; }

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 bold upper X , 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 bold upper X 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.

<<PrimarySample Public Data>>= 
int64_t lastModificationIteration = 0;

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 Method Definitions>>+=  

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

<<Enlarge MLTSampler::X if necessary and get current bold upper X Subscript i >>= 
if (index >= X.size()) X.resize(index + 1); PrimarySample &Xi = X[index];

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.

<<Reset bold upper X Subscript i if a large step took place in the meantime>>= 

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 bold upper X Subscript i ’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.

<<Apply remaining sequence of mutations to sample>>= 
Xi.Backup(); if (largeStep) { Xi.value = rng.UniformFloat(); } else { int64_t nSmall = currentIteration - Xi.lastModificationIteration; <<Apply nSmall small step mutations>> 
<<Sample the standard normal distribution upper N left-parenthesis 0 comma 1 right-parenthesis >> 
Float normalSample = Sqrt2 * ErfInv(2 * rng.UniformFloat() - 1);
<<Compute the effective standard deviation and apply perturbation to bold upper X Subscript i >> 
Float effSigma = sigma * std::sqrt((Float)nSmall); Xi.value += normalSample * effSigma; Xi.value -= std::floor(Xi.value);
} Xi.lastModificationIteration = currentIteration;

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

bold upper X prime Subscript i Baseline tilde upper N left-parenthesis bold upper X Subscript i Baseline comma sigma squared right-parenthesis comma

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

<<Apply nSmall small step mutations>>= 
<<Sample the standard normal distribution upper N left-parenthesis 0 comma 1 right-parenthesis >> 
Float normalSample = Sqrt2 * ErfInv(2 * rng.UniformFloat() - 1);
<<Compute the effective standard deviation and apply perturbation to bold upper X Subscript i >> 
Float effSigma = sigma * std::sqrt((Float)nSmall); Xi.value += normalSample * effSigma; Xi.value -= std::floor(Xi.value);

Recall the rule that when two samples from normal distributions upper N left-parenthesis mu 1 comma sigma 1 squared right-parenthesis and upper N left-parenthesis mu 2 comma sigma 2 squared right-parenthesis are added, the sum is also normally distributed with parameters upper N left-parenthesis mu 1 plus mu 2 comma sigma 1 squared plus sigma 2 squared right-parenthesis . Thus, when n perturbations are to be applied to bold upper X Subscript i , instead of performing n perturbations in sequence, it is equivalent and more efficient to directly sample

bold upper X prime Subscript i Baseline tilde upper N left-parenthesis bold upper X Subscript i Baseline comma n sigma squared right-parenthesis comma

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

p left-parenthesis x right-parenthesis equals StartFraction 1 Over StartRoot 2 pi EndRoot EndFraction normal e Superscript minus x squared slash 2

gives the following sampling method for a uniform sample xi element-of left-bracket 0 comma 1 right-parenthesis :

upper P Superscript negative 1 Baseline left-parenthesis xi right-parenthesis equals StartRoot 2 EndRoot normal e normal r normal f Superscript negative 1 Baseline left-parenthesis 2 xi minus 1 right-parenthesis comma

where normal e normal r normal f is the error function, normal e normal r normal f left-parenthesis x right-parenthesis equals 2 slash StartRoot pi EndRoot integral Subscript 0 Superscript x Baseline normal e Superscript minus x Super Superscript prime 2 Baseline normal d x Superscript prime Baseline , and normal e normal r normal f Superscript negative 1 is its inverse. The ErfInv() function, not included here, approximates normal e normal r normal f Superscript negative 1 with a polynomial.

<<Sample the standard normal distribution upper N left-parenthesis 0 comma 1 right-parenthesis >>= 
Float normalSample = Sqrt2 * ErfInv(2 * rng.UniformFloat() - 1);

We scale the resulting sample by the effective variance from Equation (16.23) and use it to perturb the sample bold upper X Subscript i before keeping only its fractional component, so that it remains in left-bracket 0 comma 1 right-parenthesis .

<<Compute the effective standard deviation and apply perturbation to bold upper X Subscript i >>= 
Float effSigma = sigma * std::sqrt((Float)nSmall); Xi.value += normalSample * effSigma; Xi.value -= std::floor(Xi.value);

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.

<<MLTSampler Method Definitions>>+=  
void MLTSampler::Reject() { for (auto &Xi : X) if (Xi.lastModificationIteration == currentIteration) Xi.Restore(); --currentIteration; }

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.

<<PrimarySample Public Methods>>= 

<<PrimarySample Public Data>>+= 
Float valueBackup = 0; int64_t modifyBackup = 0;

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 bold upper X and vertices on the camera or light subpath. With the process described above, the components bold upper X 0 comma ellipsis comma bold upper X Subscript n Baseline determine the camera subpath (for some n element-of bold upper Z Superscript ), and the remaining values bold upper X Subscript n plus 1 Baseline comma ellipsis comma bold upper X Subscript m Baseline 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 bold upper X 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.)

<<MLTSampler Method Definitions>>+= 
void MLTSampler::StartStream(int index) { streamIndex = index; sampleIndex = 0; }

<<MLTSampler Private Data>>+= 
int streamIndex, sampleIndex;

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 bold upper X are respectively used for the first dimension of each of the streams, and so forth.

<<MLTSampler Public Methods>>+= 
int GetNextIndex() { return streamIndex + streamCount * sampleIndex++; }

16.4.5 MLT Integrator

Given all of this infrastructure—an explicit representation of an n -dimensional sample  bold upper X , 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.

<<MLT Declarations>>= 
class MLTIntegrator : public Integrator { public: <<MLTIntegrator Public Methods>> 
MLTIntegrator(std::shared_ptr<const Camera> camera, int maxDepth, int nBootstrap, int nChains, int mutationsPerPixel, Float sigma, Float largeStepProbability) : camera(camera), maxDepth(maxDepth), nBootstrap(nBootstrap), nChains(nChains), mutationsPerPixel(mutationsPerPixel), sigma(sigma), largeStepProbability(largeStepProbability) { } void Render(const Scene &scene); Spectrum L(const Scene &scene, MemoryArena &arena, const std::unique_ptr<Distribution1D> &lightDistr, MLTSampler &sampler, int k, Point2f *pRaster);
private: <<MLTIntegrator Private Data>> 
std::shared_ptr<const Camera> camera; const int maxDepth; const int nBootstrap; const int mutationsPerPixel; const Float sigma, largeStepProbability; const int nChains;
};

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 upper L left-parenthesis bold upper X right-parenthesis for a vector of sample values bold upper X 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.

<<MLT Method Definitions>>= 
Spectrum MLTIntegrator::L(const Scene &scene, MemoryArena &arena, const std::unique_ptr<Distribution1D> &lightDistr, MLTSampler &sampler, int depth, Point2f *pRaster) { sampler.StartStream(cameraStreamIndex); <<Determine the number of available strategies and pick a specific one>> 
int s, t, nStrategies; if (depth == 0) { nStrategies = 1; s = 0; t = 2; } else { nStrategies = depth + 2; s = std::min((int)(sampler.Get1D() * nStrategies), nStrategies - 1); t = nStrategies - s; }
<<Generate a camera subpath with exactly t vertices>> 
Vertex *cameraVertices = arena.Alloc<Vertex>(t); Bounds2f sampleBounds = (Bounds2f)camera->film->GetSampleBounds(); *pRaster = sampleBounds.Lerp(sampler.Get2D()); if (GenerateCameraSubpath(scene, sampler, arena, t, *camera, *pRaster, cameraVertices) != t) return Spectrum(0.f);
<<Generate a light subpath with exactly s vertices>> 
sampler.StartStream(lightStreamIndex); Vertex *lightVertices = arena.Alloc<Vertex>(s); if (GenerateLightSubpath(scene, sampler, arena, s, cameraVertices[0].time(), *lightDistr, lightVertices) != s) return Spectrum(0.f);
<<Execute connection strategy and return the radiance estimate>> 
sampler.StartStream(connectionStreamIndex); return ConnectBDPT(scene, lightVertices, cameraVertices, s, t, *lightDistr, *camera, sampler, pRaster) * nStrategies;
}

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 s equals 1 or t equals 1 (refer to Section 16.3.3 for details).

<<MLTSampler Constants>>= 
static const int cameraStreamIndex = 0; static const int lightStreamIndex = 1; static const int connectionStreamIndex = 2; static const int nSampleStreams = 3;

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 monospace d monospace e monospace p monospace t monospace h plus 2 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.

<<Determine the number of available strategies and pick a specific one>>= 
int s, t, nStrategies; if (depth == 0) { nStrategies = 1; s = 0; t = 2; } else { nStrategies = depth + 2; s = std::min((int)(sampler.Get1D() * nStrategies), nStrategies - 1); t = nStrategies - s; }

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 left-bracket 0 comma 1 right-parenthesis squared , maps it to raster coordinates, and tries to generate a corresponding camera subpath with exactly t vertices, failing with a 0-valued estimate for upper L when this was not possible.

<<Generate a camera subpath with exactly t vertices>>= 
Vertex *cameraVertices = arena.Alloc<Vertex>(t); Bounds2f sampleBounds = (Bounds2f)camera->film->GetSampleBounds(); *pRaster = sampleBounds.Lerp(sampler.Get2D()); if (GenerateCameraSubpath(scene, sampler, arena, t, *camera, *pRaster, cameraVertices) != t) return Spectrum(0.f);

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

<<MLTIntegrator Private Data>>= 
std::shared_ptr<const Camera> camera;

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.

<<Generate a light subpath with exactly s vertices>>= 
sampler.StartStream(lightStreamIndex); Vertex *lightVertices = arena.Alloc<Vertex>(s); if (GenerateLightSubpath(scene, sampler, arena, s, cameraVertices[0].time(), *lightDistr, lightVertices) != s) return Spectrum(0.f);

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

<<Execute connection strategy and return the radiance estimate>>= 
sampler.StartStream(connectionStreamIndex); return ConnectBDPT(scene, lightVertices, cameraVertices, s, t, *lightDistr, *camera, sampler, pRaster) * 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 b equals integral upper I left-parenthesis bold upper X right-parenthesis normal d normal upper Omega , 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.

<<MLT Method Definitions>>+= 
void MLTIntegrator::Render(const Scene &scene) { std::unique_ptr<Distribution1D> lightDistr = ComputeLightPowerDistribution(scene); <<Generate bootstrap samples and compute normalization constant b >> 
int nBootstrapSamples = nBootstrap * (maxDepth + 1); std::vector<Float> bootstrapWeights(nBootstrapSamples, 0); std::vector<MemoryArena> bootstrapThreadArenas(MaxThreadIndex()); ParallelFor( [&](int i) { <<Generate ith bootstrap sample>> 
MemoryArena &arena = bootstrapThreadArenas[ThreadIndex]; for (int depth = 0; depth <= maxDepth; ++depth) { int rngIndex = i * (maxDepth + 1) + depth; MLTSampler sampler(mutationsPerPixel, rngIndex, sigma, largeStepProbability, nSampleStreams); Point2f pRaster; bootstrapWeights[rngIndex] = L(scene, arena, lightDistr, sampler, depth, &pRaster).y(); arena.Reset(); }
}, nBootstrap, 4096); Distribution1D bootstrap(&bootstrapWeights[0], nBootstrapSamples); Float b = bootstrap.funcInt * (maxDepth + 1);
<<Run nChains Markov chains in parallel>> 
Film &film = *camera->film; int64_t nTotalMutations = (int64_t)mutationsPerPixel * (int64_t)film.GetSampleBounds().Area(); ParallelFor( [&](int i) { int64_t nChainMutations = std::min((i + 1) * nTotalMutations / nChains, nTotalMutations) - i * nTotalMutations / nChains; <<Follow ith Markov chain for nChainMutations>> 
MemoryArena arena; <<Select initial state from the set of bootstrap samples>> 
RNG rng(i); int bootstrapIndex = bootstrap.SampleDiscrete(rng.UniformFloat()); int depth = bootstrapIndex % (maxDepth + 1);
<<Initialize local variables for selected state>> 
MLTSampler sampler(mutationsPerPixel, bootstrapIndex, sigma, largeStepProbability, nSampleStreams); Point2f pCurrent; Spectrum LCurrent = L(scene, arena, lightDistr, sampler, depth, &pCurrent);
<<Run the Markov chain for nChainMutations steps>> 
for (int64_t j = 0; j < nChainMutations; ++j) { sampler.StartIteration(); Point2f pProposed; Spectrum LProposed = L(scene, arena, lightDistr, sampler, depth, &pProposed); <<Compute acceptance probability for proposed sample>> 
Float accept = std::min((Float)1, LProposed.y() / LCurrent.y());
<<Splat both current and proposed samples to film>> 
if (accept > 0) film.AddSplat(pProposed, LProposed * accept / LProposed.y()); film.AddSplat(pCurrent, LCurrent * (1 - accept) / LCurrent.y());
<<Accept or reject the proposal>> 
if (rng.UniformFloat() < accept) { pCurrent = pProposed; LCurrent = LProposed; sampler.Accept(); } else sampler.Reject();
}
}, nChains);
<<Store final image computed with MLT>>  }

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 monospace m monospace a monospace x monospace upper D monospace e monospace p monospace t monospace h monospace plus monospace 1 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.

<<Generate bootstrap samples and compute normalization constant b >>= 
int nBootstrapSamples = nBootstrap * (maxDepth + 1); std::vector<Float> bootstrapWeights(nBootstrapSamples, 0); std::vector<MemoryArena> bootstrapThreadArenas(MaxThreadIndex()); ParallelFor( [&](int i) { <<Generate ith bootstrap sample>> 
MemoryArena &arena = bootstrapThreadArenas[ThreadIndex]; for (int depth = 0; depth <= maxDepth; ++depth) { int rngIndex = i * (maxDepth + 1) + depth; MLTSampler sampler(mutationsPerPixel, rngIndex, sigma, largeStepProbability, nSampleStreams); Point2f pRaster; bootstrapWeights[rngIndex] = L(scene, arena, lightDistr, sampler, depth, &pRaster).y(); arena.Reset(); }
}, nBootstrap, 4096); Distribution1D bootstrap(&bootstrapWeights[0], nBootstrapSamples); Float b = bootstrap.funcInt * (maxDepth + 1);

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

<<MLTIntegrator Private Data>>+=  
const int maxDepth; const int nBootstrap;

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

<<Generate ith bootstrap sample>>= 
MemoryArena &arena = bootstrapThreadArenas[ThreadIndex]; for (int depth = 0; depth <= maxDepth; ++depth) { int rngIndex = i * (maxDepth + 1) + depth; MLTSampler sampler(mutationsPerPixel, rngIndex, sigma, largeStepProbability, nSampleStreams); Point2f pRaster; bootstrapWeights[rngIndex] = L(scene, arena, lightDistr, sampler, depth, &pRaster).y(); arena.Reset(); }

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.

<<MLTIntegrator Private Data>>+=  
const int mutationsPerPixel; const Float sigma, largeStepProbability;

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.

<<Run nChains Markov chains in parallel>>= 
Film &film = *camera->film; int64_t nTotalMutations = (int64_t)mutationsPerPixel * (int64_t)film.GetSampleBounds().Area(); ParallelFor( [&](int i) { int64_t nChainMutations = std::min((i + 1) * nTotalMutations / nChains, nTotalMutations) - i * nTotalMutations / nChains; <<Follow ith Markov chain for nChainMutations>> 
MemoryArena arena; <<Select initial state from the set of bootstrap samples>> 
RNG rng(i); int bootstrapIndex = bootstrap.SampleDiscrete(rng.UniformFloat()); int depth = bootstrapIndex % (maxDepth + 1);
<<Initialize local variables for selected state>> 
MLTSampler sampler(mutationsPerPixel, bootstrapIndex, sigma, largeStepProbability, nSampleStreams); Point2f pCurrent; Spectrum LCurrent = L(scene, arena, lightDistr, sampler, depth, &pCurrent);
<<Run the Markov chain for nChainMutations steps>> 
for (int64_t j = 0; j < nChainMutations; ++j) { sampler.StartIteration(); Point2f pProposed; Spectrum LProposed = L(scene, arena, lightDistr, sampler, depth, &pProposed); <<Compute acceptance probability for proposed sample>> 
Float accept = std::min((Float)1, LProposed.y() / LCurrent.y());
<<Splat both current and proposed samples to film>> 
if (accept > 0) film.AddSplat(pProposed, LProposed * accept / LProposed.y()); film.AddSplat(pCurrent, LCurrent * (1 - accept) / LCurrent.y());
<<Accept or reject the proposal>> 
if (rng.UniformFloat() < accept) { pCurrent = pProposed; LCurrent = LProposed; sampler.Accept(); } else sampler.Reject();
}
}, nChains);

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 amount of time.

<<MLTIntegrator Private Data>>+= 
const int nChains;

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.

<<Follow ith Markov chain for nChainMutations>>= 
MemoryArena arena; <<Select initial state from the set of bootstrap samples>> 
RNG rng(i); int bootstrapIndex = bootstrap.SampleDiscrete(rng.UniformFloat()); int depth = bootstrapIndex % (maxDepth + 1);
<<Initialize local variables for selected state>> 
MLTSampler sampler(mutationsPerPixel, bootstrapIndex, sigma, largeStepProbability, nSampleStreams); Point2f pCurrent; Spectrum LCurrent = L(scene, arena, lightDistr, sampler, depth, &pCurrent);
<<Run the Markov chain for nChainMutations steps>> 
for (int64_t j = 0; j < nChainMutations; ++j) { sampler.StartIteration(); Point2f pProposed; Spectrum LProposed = L(scene, arena, lightDistr, sampler, depth, &pProposed); <<Compute acceptance probability for proposed sample>> 
Float accept = std::min((Float)1, LProposed.y() / LCurrent.y());
<<Splat both current and proposed samples to film>> 
if (accept > 0) film.AddSplat(pProposed, LProposed * accept / LProposed.y()); film.AddSplat(pCurrent, LCurrent * (1 - accept) / LCurrent.y());
<<Accept or reject the proposal>> 
if (rng.UniformFloat() < accept) { pCurrent = pProposed; LCurrent = LProposed; sampler.Accept(); } else sampler.Reject();
}

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.

<<Select initial state from the set of bootstrap samples>>= 
RNG rng(i); int bootstrapIndex = bootstrap.SampleDiscrete(rng.UniformFloat()); int depth = bootstrapIndex % (maxDepth + 1);

Having chosen the bootstrap sample, we must now obtain the corresponding primary sample space vector bold upper X . 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 exact same sampler that originally produced the sampled entry of the bootstrap distribution.

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

<<Initialize local variables for selected state>>= 
MLTSampler sampler(mutationsPerPixel, bootstrapIndex, sigma, largeStepProbability, nSampleStreams); Point2f pCurrent; Spectrum LCurrent = L(scene, arena, lightDistr, sampler, depth, &pCurrent);

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.

<<Run the Markov chain for nChainMutations steps>>= 
for (int64_t j = 0; j < nChainMutations; ++j) { sampler.StartIteration(); Point2f pProposed; Spectrum LProposed = L(scene, arena, lightDistr, sampler, depth, &pProposed); <<Compute acceptance probability for proposed sample>> 
Float accept = std::min((Float)1, LProposed.y() / LCurrent.y());
<<Splat both current and proposed samples to film>> 
if (accept > 0) film.AddSplat(pProposed, LProposed * accept / LProposed.y()); film.AddSplat(pCurrent, LCurrent * (1 - accept) / LCurrent.y());
<<Accept or reject the proposal>> 
if (rng.UniformFloat() < accept) { pCurrent = pProposed; LCurrent = LProposed; sampler.Accept(); } else sampler.Reject();
}

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 upper L left-parenthesis bold upper X right-parenthesis 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.

<<Compute acceptance probability for proposed sample>>= 
Float accept = std::min((Float)1, LProposed.y() / LCurrent.y());

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.

<<Splat both current and proposed samples to film>>= 
if (accept > 0) film.AddSplat(pProposed, LProposed * accept / LProposed.y()); film.AddSplat(pCurrent, LCurrent * (1 - accept) / LCurrent.y());

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.

<<Accept or reject the proposal>>= 
if (rng.UniformFloat() < accept) { pCurrent = pProposed; LCurrent = LProposed; sampler.Accept(); } else sampler.Reject();

Metropolis sampling only considers the relative frequency of samples and cannot create an image that is correctly scaled in absolute terms; hence the value b 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 b when writing the image to convert to actual incident radiance on the film.

<<Store final image computed with MLT>>=