An embedded domain specific sampling language for Monte Carlo rendering

Getting Started

Language Overview


Aether is a domain specific sampling language. It is embedded in C++ and was designed primarily for use in Monte Carlo rendering applications. The user writes sampling code, then at compile time Aether automatically derives the corresponding PDF code. The generated PDF code is consistent with the provided sampling code, making it easier to write correct Monte Carlo estimators.

To learn more, see the language overview or tutorials; read the paper; or check out the code for Aether or our Mitsuba/Aether renderer, Yotsuba.

Getting Started

Downloading Aether

Aether is available here.

Compiling Aether

Aether is header only but uses features of C++11/14/17 that require a recent version of clang to compile. It has been tested with clang 3.7. The language makes extensive use of template metaprogramming and requires the following settings when compiling:



All Aether code is wrapped in the aether namespace.

Language Overview

An Example: Estimating Irradiance at a Point

Suppose we want to compute the irradiance at a fixed point by using Monte Carlo integration to evaluate a hemispherical integral.

Sampling the Hemisphere

We want to importance sample from a cosine weighted hemisphere:

// Declare symbolic uniform random variables
variable<1> u1;
variable<2> u2;

auto r = sqrt(u1);
auto phi = 2_l * pi * u2; // 2_l is a literal for 2.0
auto cosHemisphere = random_vector<2>(
  r * cos(phi),
  r * sin(phi),
  sqrt(1 - u1)

The above code looks similar to regular numeric sampling code, but all of the expressions are symbolic. It starts with symbolic uniform random variables u1 and u2 and defines the symbolic expression to transform them into a random variable (cosHemisphere) representing directions on the hemisphere. The random_vector<N>() function constructs a vector of random variables, where 2 is the number of uniform random variables on which it depends. In this example, there are 2 uniforms — u1 and u2 — and the random vector has 3 outputs, because directions are two-dimensional but encoded with 3 coordinates.

Generating Samples

Since cosHemisphere is symbolic, at this stage, nothing has actually been sampled. Every random vector in the language has a Sample method. If we want to generate samples from cosHemisphere, we call Sample with 2 numbers between 0 and 1, one each for u1 and u2:

auto dir = cosHemisphere.Sample(0.25f, 0.75f);

It is important to note that the numbers provided to Sample() can be generated however we choose. Typically they will be generated by some pseudo-random number generator or quasi-Monte Carlo sequence, for example:

auto uv = someRNG.Next2D();
auto dir = cosHemisphere.Sample(uv[0], uv[1]);

Computing the PDF

At compile time, Aether automatically computes the symbolic Jacobian determinant and symbolic inverse, both of which are used when computing the PDF of samples. cosHemisphere has an automatically provided Pdf method.

auto dirPdf = cosHemisphere.Pdf(dir);

This method can compute not only the PDF of directions sampled by cosHemisphere but any given direction. This becomes particularly valuable when combining samples.

Computing an Estimate

Using Aether’s Pdf() method, computing an estimate is simple:

// User provided uniform sampler
MyRng rng;

Spectrum total(0);
int N = 10000;

auto myIntegrand = [](const auto& sample) -> Spectrum {
  // user provided regular C++ code to compute Li * |cos(theta)|

for (int i = 0; i < N; i++) {
  // Draw a sample
  auto xCos = cosHemisphere.Sample(rng(), rng());

  // Compute the integrand
  auto f = myIntegrand(xCos);

  // Compute the PDF
  auto p = cosHemisphere.Pdf(xCos);

  // Add f(xCos) / p(xCos) to running total
  total += f / p;

// Compute the final estimate
auto estimate = total / N;

This estimator is correct by construction because Aether uses symbolic differentiation to ensure that the sampling and PDF code are consistent.

Discrete Random Variables

Discrete random variables in Aether are constructed from standard C++ containers. They support the same Sample() and Pdf() methods as the continuous random variables.

std::vector<Emitter*> emitters = scene.getEmitters();
// Create a uniform discrete distribution
auto emitterDiscrete = discrete(emitters);

// Sample an emitter
auto em = emitterDiscrete.Sample(rng());

// Compute its probability
auto p = emitterDiscrete.Pdf(em);


Discrete random variables in Aether are constructed from standard C++ containers. They support the same Sample() and Pdf() methods as the continuous random variables.

std::vector<Emitter*> emitters = scene.getEmitters();
// Create a uniform discrete distribution
auto emitterDiscrete = discrete(emitters);

// Sample an emitter
auto em = emitterDiscrete.Sample(rng());

// Compute its probability
auto p = emitterDiscrete.Pdf(em);

Multiple Importance Sampling

Multiple importance sampling requires evaluating several probability densities for each sample drawn. Suppose we wish to use MIS to combine the cosine hemisphere samples with samples from an analogous uniform hemisphere random variable uniformHemisphere. A key feature of Aether that enables this is that the Pdf() method of random variables accepts as argument not just samples generated by its own Sample() method, but any sample.

As a result, evaluating all four densities is easy. Denoting the cosine weighted sample as xCos and the uniform sample as xUniform, we simply compute:

// PDF of xCos if sampled by cosHemisphere
auto pCos = xCos.Pdf(xCos);
// PDF of xCos if sampled by uniformHemisphere
auto pUniformCos = xUniform.Pdf(xCos);
// PDF of xUniform if sampled by uniformHemisphere
auto pUniform = xUniform.Pdf(xUniform);
// PDF of xUniform if sampled by cosHemisphere
auto pCosUniform = xCos.Pdf(xUniform);

All samples in Aether also store the random variable from which they were sampled, so we can use the Pdf() method on xCos and xUniform directly to compute the PDFs. However, no extra storage is needed at run time because all such dependencies are resolved statically.

If the sample provided to Pdf() is not of the correct dimension or domain, the inverse will return invalid uniforms (outside [0, 1]) and Pdf() simply returns 0. The programmer does not need to manually check for these cases.


To simplify the process of combining potentially arbitrary numbers of samples, Aether provides the combine() function, which accepts a user defined combining heuristic and any number of samples. The above example can be written more concisely as:

Estimator<Spectrum> myEstimator(myIntegrand);

// User provided code to compute the MIS weight
struct PowerHeuristic {
  auto operator()(float pdfA, float pdfB) const {
    // e.g. the power heuristic:
    return (pdfA * pdfA) / (pdfA * pdfA + pdfB * pdfB);

// Combine samples with power heuristic
auto combined = combine<PowerHeuristic>(xCos, xUniform);

// Accumulate weighted integrand values 
total += myEstimator(combined);

Interfacing with Deterministic Code

Aether is designed to be usable with external code such as ray casting engines. This creates the need to compute densities that depend on data computed outside the language, for which symbolic descriptions are unavailable. To balance the need for such an interface with the need for symbolic derivation, we require that the data coming from the deterministic external code be constant in the neighborhood of a sample. This means, for example, that ray casting cannot directly return an intersection point, but should instead return the triangle’s vertices, which are constant in a neighborhood of the intersection, and that the intersection point coordinates must be computed in our language. This ensures that proper derivatives, inverses and densities can be derived symbolically.

Aether provides the ConstantCall mechanism for calling external deterministic code. A triangle intersection is implemented as

// Intersect ray (p, dir) with the scene and expect a constant as a result
Intersection its = context.ConstantCall(raycaster, Ray(p.Value(), dir.Value()));

// Compute ray-triangle intersection in Aether
auto v0 = constant(its.v0);
auto v1 = constant(its.v1);
auto v2 = constant(its.v2);
auto e1 = v0 - v2;
auto e2 = v1 - v2;
auto N = cross(e1, e2); 
auto t = dot(v0 - p, N) / dot(dir, N);
return  p + t * dir;

The Intersection object returned by the ray casting engine includes the vertices (and other relevant information) of the intersected object. We first cast each vertex to a constant (constant(its.v0), etc.) and implement a standard ray-triangle intersection in Aether to obtain the final intersection point.

Random Sequences

The random vector introduced above has a fixed size and is designed for situations when all its coordinates are sampled at once. Random sequences instead support the creation of incremental sequences of generic random variables. Each element of a random sequence is sampled from a strategy. We use this type to represent light paths. The type of data stored in the sequence is user provided, e.g., a Vertex type representing a point on a surface.

Consider sampling a 2-vertex path segment starting from an emitter. This is implemented by a random sequence of two strategies, SamplePointOnLightStrategy (defined above), and SampleHemiAndIntersectStrategy, a strategy that picks a cosine weighted direction, traces a ray starting at the previous path vertex, and computes the intersection:

// Strategy for sampling a cosine-weighted direction
// and intersecting it with the scene
struct SampleHemiAndIntersectStrategy {
  template <typename T>
  auto operator()(Context<T>& context,
                  const RandomSequence<Vertex>& path,
                  MyRng& rng,
                  Raycaster& raycaster) const {
    // Sample uniforms between 0 and 1
    auto uv = context.Uniform2D(rng);

    // Sample the outgoing direction
    // cosHemisphere is defined as above
    auto dir = cosHemisphere.Sample(uv[0], uv[1]);

    // Get the previous Vertex from RandomSequence
    auto p = path.Back();
    // Compute intersection for ray (p, dir) as above
    // ...

// Create an initially empty path
RandomSequence<Vertex> path;

// Append a point on a light
SamplePointOnLightStrategy samplePointOnLight;
path.Append(samplePointOnLight, rng, emitters);

// Append an intersection point
SampleHemiAndIntersectStrategy sampleHemiAndIntersect;
path.Append(sampleHemiAndIntersect, rng, raycaster);

// Actually sample the strategies

The fundamental operation of a random sequence is to Append new elements to the sequence. When path.Append(samplePointOnLight, ...) is called, the samplePointOnLight strategy is stored, without being sampled. Random sequences are evaluated lazily. When Sample() is called, each strategy is then sampled in sequence with path as an argument (along with any other provided arguments). The result is a new Vertex sample, which is stored along with the strategy from which it was sampled.

Like the other random variables in Aether, random sequences offer a Pdf() method. Like all other random variables, the random sequence Pdf() method can be used to compute the PDF of any sequence, not just itself, and hence random sequences can be combined with MIS just as easily.

Aether also provides the functions slice, concat, and reverse for extracting subsequences, concatenating sequences, and reversing their order.

The separation of strategies and random sequences has the additional benefit that strategies can be built to be orthogonal and easily reused: they are not tied to a specific random sequence nor to a particular algorithm. For instance, in the above example, changing the cosine weighted hemisphere sampling to uniform sampling would be as simple as defining a new strategy based on uniformHemisphere and appending that instead.

Rendering Algorithm Tutorials

Yotsuba is a Mitsuba based renderer that uses Aether to generate samples and compute PDFs. It is available here and contains implementations of many standard rendering algorithms.

The following sections describe how to write a rendering algorithm in Yotsuba.

Writing a Direct Lighting Integrator in Yotsuba

We first need to create a new integrator. See Mitsuba’s documentation for more details on how to do this. The main logic of our direct lighting integrator will be implemented in the Li() method:

Spectrum Li(const RayDifferential &ray, RadianceQueryRecord &rRec) const {
  // Our direct lighting integrator logic

Writing rendering algorithms with Aether is different from traditional approaches like in pbrt and Mitsuba.

Create a Path

Aether provides a path data structure (RandomSequence) that is used as the basic building block of all rendering algorithms. We sample an entire light path by repeatedly appending new sampled vertices to this data structure. But to start, we simply create an empty path:

// Create camera path
RandomSequence<Vertex> camPath;

Append the Camera Position

In order to append new vertices to our path, All vertices in a RandomSequence are sampled from strategies. In order to append the camera position as the first vertex in our path, we need to create a strategy:

struct sample_position_on_sensor_t {
  template <typename T>
  auto operator()(Context<T>& context, const RandomSequence<Vertex>& seq, const Ray& ray) const {
    return sample(constant(to_vector3(ray.o)));

All strategies require a templated Context<T>& context object as first argument. This is necessary for computing the PDF of strategies, which may contain both discrete and continuous distributions. More details on this are available in the paper.

The second argument is the RandomSequence to which this strategy will be appended. Additional arguments can be provided as needed. In this case, we provide the initial camera provided by Mitsuba.

The body of the strategy defines a random variable for sampling the next vertex in the path. In this case, we want to simply sample the camera vertex, which is given by the ray’s origin (ray.o). to_vector3 converts the origin from Mitsuba’s data type to Aether’s. constant creates a constant random variable, with the value of the given point. And the result if then provided to sample as the result of the function.

It is important to note that sampling the camera position is a special case in Aether. The camera position (and the initial sampled ray) is not actually sampled randomly, it has a fixed position. Formally, it is defined by a Dirac delta function. Aether does not currently provide a construct for handling Diracs. We assume the ratio between position density and the integrand is 1. A constant random variable in Aether has a PDF of 1, which correctly accounts for the ratio between density and integrand, but is not formally a Dirac. So the use of constant here is slightly misleading, but is only the case when sampling these camera paths.

We first create an instance of the strategy:

Node<sample_position_on_sensor_t> sampleCamPos;

Then we append the strategy to our path:

// Append the camera position
camPath.Append(sampleCamPos, ray);

The Append method simply stores the provided strategy along with its arguments. When the vertices are actually sampled, Aether will create a context object and call sampleCamPos(context, camPath, ray); to actually generate the new vertex.

Append the Primary Intersection

At this point, camPath contains 1 strategy for the camera position. We now extend it to obtain the primary intersection vertex. As before, we create a strategy:

struct sample_primary_intersection_t {
  Raycaster& raycaster;

  template <typename T>
  auto operator()(Context<T>& context, const RandomSequence<Vertex>& seq, const Ray& ray) const {
    const Intersection its = context.constant_call(raycaster, ray);
    auto pt = constant(its.p);

    return sample(
      , intersection_ = its

In this strategy, we need to cast a ray into the scene. We wrap the call to raycaster in context.constant_call(). This tells Aether that we are calling deterministic code outside the language. The raycaster object is provided by the user.

The line, intersection_ = its, tells Aether that we want to store the intersection information for this vertex. This information will be necessary when we compute the integrand.

As before, we create an instance of the strategy, passing a raycaster object to the constructor:

Node<sample_primary_intersection_t> samplePrimaryIntersection{raycaster};

Then we append the strategy to our path:

// Append the primary intersection
camPath.Append(samplePrimaryIntersection, ray);

Append a Light Vertex

To complete our direct lighting path, we want to sample a vertex on a discretely chosen light source.

struct sample_position_on_emitter_t {
  template <typename T>
  auto operator()(Context<T>& context, const RandomSequence<Vertex>& seq, UniDist& uniDist, const ref_vector<Emitter> &emitters) const {
		auto emitters_rv = discrete_dynamic(emitters);
		auto emitter = context.Sample(emitters_rv, context.Uniform1D(uniDist));
		auto sampler = emitter->makeSampler();
		auto pt = sampler.Sample(context, uniDist);

		return sample(
			, emitter_ = emitter.get()

This strategy introduces several new constructs. discrete_dynamic creates a new uniform discrete random variable over the provided emitters. To actually sample an emitter from this random variable, we call context.Sample(emitters_rv, context.Uniform1D(uniDist)). For all discrete random variables inside strategies, we use context.Sample() instead of calling emitters_rv.Sample() directly. This is required for computing the PDF over the strategy and otherwise does not impact the user. For the same reason, we call context.Uniform1D(uniDist)) to obtain a single uniform between 0 and 1, instead of calling uniDist directly. Next, we call emitter->makeSampler(), which creates a random variable for actually sampling the new position on the light. The actual definition of makeSampler is given below. We then actually generate the sample vertex by calling sampler.Sample() and return it along with the emitter itself, which will be used when computing the integrand.

As before, we create an instance of the strategy:

Node<sample_position_on_emitter_t> sampleEmtPosition;

Then we append the strategy to our path, passing a uniform random number generator (uniDist) and the scene’s list of emitters as arguments:

// Append an emitter vertex 
camPath.Append(sampleEmtPosition, uniDist, scene->emitters);

Defining Different Emitter Samplers

When we sample emitter, we don’t know what type of emitter it is. It might be a triangle mesh emitter, a sphere emitter, an environment map, or something else. In typical renderers, each emitter type is derived from some base class and overrides a virtual method for sampling a point on its surface. In Aether, we require that the result when sampling the point is a random variable in the language. All random variables are defined by a compile time symbolic expression (a template expression), which is encoded in its type. As such, random variables defined by different expressions have different types. So the method that returns the random variable must have a different return type for each of them. This is an issue because virtual methods require a consistent return type.

In order to support different return types, Aether provides a composite random variable construct. We create a compile time composite random variable over all the possible return values, then the virtual function selects at run time which of the return values should be selected. The composite random variable acts as a kind of table, where the possible values are supplied at compile time and the actual value is then selected at run time.

Suppose we have a single emitter type, a triangle emitter. We first need to define its sample method:

struct trimesh_emitter {
  const std::vector<const Triangle*> *trianglePtrs;
  const std::vector<Real> *triangleWeights;
  const Point *positions;
  const Normal *normals;

  template <typename ContextType>
  auto SamplePosition(Context<ContextType>& context, UniDist& uniDist) const {
    // Create a discrete random variable over all triangles
    auto rv = discrete_dynamic(*trianglePtrs, *triangleWeights);

    // Pick 1 of the triangles
    auto trianglePtr = context.Sample(rv, context.Uniform1D(uniDist));

    const Vector3 v0 = to_vector3(positions[trianglePtr->idx[0]]);
    const Vector3 v1 = to_vector3(positions[trianglePtr->idx[1]]);
    const Vector3 v2 = to_vector3(positions[trianglePtr->idx[2]]);

    auto uv = context.Uniform2D(uniDist);
    auto pt = uniform_triangle(v0, v1, v2).Sample(uv[0], uv[1]);

    // Fill intersection information for integrand evaluation
    Intersection its;
    its.p = to_point(pt.Value());
    its.geoFrame = Frame(to_normal((v1 - v0).cross(v2 - v0).normalized()));
    // Phong normal interpolation
    const Normal n0 = normals[trianglePtr->idx[0]];
    const Normal n1 = normals[trianglePtr->idx[1]];
    const Normal n2 = normals[trianglePtr->idx[2]];
    const Float a = sqrt(uv[0]);
    const Float b0 = 1.f - a;
    const Float b1 = uv[1] * a;
    its.shFrame = Frame(b0 * n0 + b1 * n1 + (1.f - b0 - b1) * n2);

    return sample(
      , v0_ = v0
      , v1_ = v1
      , v2_ = v2
      , intersection_ = its

The SamplePosition method first discretely picks an individual triangle from the emitter’s triangle mesh. It then simply samples a point uniformly on the triangle (the uniform_triangle function is part of the standard library of Aether) and fills in some intersection information needed for integrand evaluation.

We need to create a composite random variable with all the possible emitter types. In this case we just have one:

using emitter_composite_t = CompositeRandomVar<trimesh_emitter>;

In order for Aether to know which method to call when the composite random variable is sampled, we need the following template wrapper:

template <>
struct SampleCall<mitsuba::trimesh_emitter> {
  template <int I, typename Cond, typename ContextType>
  auto operator()(const mitsuba::trimesh_emitter& emitter, _int<I>, Cond, aether::Context<ContextType>& context, mitsuba::UniDist &uniDist) const {
    return emitter.SamplePosition(context, uniDist);

This functor just forwards the given arguments to the emitter’s actual sampling method. The other arguments to the method are internal to Aether and aren’t relevant to the user.

We need to create a virtual method on Mitsuba’s emitter to return the composite:

virtual emitter_composite_t makeSampler() const {
  return {};

Finally, we override this virtual method in Mitsuba’s area light definition to return a trimesh_emitter:

emitter_composite_t makeSampler() const {
  ref<TriMesh> triMesh = m_shape->createTriMesh();
  std::vector<Real> triangleWeights;
  std::vector<const Triangle*> trianglePtrs;
  size_t triangleCount = triMesh->getTriangleCount();
  const Triangle* trianglePtr = triMesh->getTriangles();
  const Point* vertexPositions = triMesh->getVertexPositions();
  const Normal* vertexNormals = triMesh->getVertexNormals();
  for (size_t i = 0; i < triangleCount; i++) {
  for (size_t i = 0; i < triangleCount; i++) {
    // Sample according to area
  return trimesh_emitter{&trianglePtrs, &triangleWeights, vertexPositions, vertexNormals};

This method collects the necessary data for the triangle vertices, as required to construct a trimesh_emitter. It’s more efficient to construct this trimesh_emitter once, then store it as a member of area light, and just return the member variable from makeSampler(), but we include the full definition here for simplicity.

When this method is called at run time, it returns a trimesh_emitter, which selects the appropriate table entry from the composite random variable.

Compute an Estimate

At this point, the path contains 3 strategies and represents a direct lighting path. Random sequences in Aether are sampled lazily so nothing has actually been sampled yet. To actually sample the vertices of our direct lighting path, we can simply call:


We are now ready to compute an estimate for camPath. We need to compute the integrand. Aether requires the user to provide an integrand function. It should compute the radiance along a given path:

Spectrum integrand(const RandomSequence<Vertex>& path) {
  // Compute the radiance along path

Aether provides additional methods for accessing data in a RandomSequence:

See the API for more details.

To complete the estimate, we need to compute the PDF. Aether automatically generates code for computing the PDF. The PDF of our path can be obtained by calling:


So the final estimate is simply:

auto estimate = integrand(camPath) / camPath.Pdf();

Writing a Path Tracer in Mitsuba

One of the benefits of Aether is that its constructs are designed to be reusable. Strategies can easily be reused in different algorithms so if we want to extend our direct lighting integrator and write a path tracer, we can use much of what we have already done.

We start by creating a path and appending the camera position and the primary intersection, reusing the strategies from above:

// Create camera subpath
RandomSequence<Vertex> camPath;
// Append the camera position
camPath.Append(sampleCamPos, ray);
// Append the primary intersection
camPath.Append(samplePrimaryIntersection, ray);

Append BSDF Vertices

Now we want to extend the path by sampling the BSDF and intersecting with the scene. As before, we create a strategy:

struct sample_bsdf_t {
  Raycaster& raycaster;

  template <typename T>
  auto operator()(Context<T>& context, const RandomSequence<Vertex>& seq, UniDist& uniDist) const {
    // Create coordinate basis at previous intersection point
    const Intersection &currentIts = seq.Back().Get(intersection_);
    auto shadingFrame = make_random_frame(currentIts);

    // Get direction to previous vertex
    auto directionToPreviousWorldVal = (seq[seq.Size() - 2].Value() - seq[seq.Size() - 1].Value()).normalized();
    auto directionToPreviousWorld = constant(directionToPreviousWorldVal);

    // Convert to local coordinates
    auto directionToPreviousLocal = to_local(shadingFrame, directionToPreviousWorld);

    // Sample a new outgoing direction
    const BSDF *bsdf = currentIts.getBSDF();
    auto bsdfSampler = bsdf->makeSampler(currentIts.toLocal(to_vector(directionToPreviousWorldVal)), currentIts);

    // Convert outgoing direction to world coordinates
    auto directionToNextLocal = bsdfSampler.Sample(directionToPreviousLocal, context, uniDist);
    auto directionToNextWorld = to_world(shadingFrame, directionToNextLocal);

    auto p = constant(seq.Back().Value());

    // Obtain next intersection point
    const Ray ray(to_point(p.Value()), to_vector(directionToNextWorld.Value()));
    const Intersection nextIts = context.constant_call(raycaster, ray);
    auto v0 = to_vector3(nextIts.triangle_v0);
    auto v1 = to_vector3(nextIts.triangle_v1);
    auto v2 = to_vector3(nextIts.triangle_v2);

    // Recompute intersection point in Aether
    auto pt = intersect(v0, v1, v2, p, directionToNextWorld);

    return optional_sample(
      , pt
      , v0_ = v0
      , v1_ = v1
      , v2_ = v2
      , intersection_ = nextIts
      , emitter_ = emitter

This strategy introduces two new concepts. After casting the ray into the scene, we recompute the intersection point in Aether. The reason for this is that in order for Aether to automatically compute the PDF, it needs derivatives of the intersection point with respect to the uniforms that were used to sample it. In order to compute the derivatives, the intersection point must be a symbolic random variable. But the intersection point returned by the raycaster is a numerical 3D point. To obtain a symbolic version of the intersection point, the raycaster must return the vertices (v0, v1, v2) of the intersected triangle, which are constant in the neighborhood of the intersection point, and we then recompute the ray-triangle intersection in Aether, using the provided intersect function. The vertices (v0, v1, v2) are also stored as additional information with the sample. This is also necessary for computing the PDF.

In previous strategies, we returned the result in a call to sample. Here, we instead return the result in a call to optional_sample. This is used to indicate to Aether if the sample is actually valid, i.e. that the ray actually intersected the scene. The first argument to optional_sample should be true for valid samples and false for invalid samples.

As before, we create an instance of the strategy, passing a raycaster object to the constructor:

Node<sample_bsdf_t> sampleBSDF{raycaster};

We can then append the strategy as normal, but here we want to sample the BSDF multiple times to extend the path:

// Extend by successive BSDF sampling until max depth
for(; camPath.Size() <= maxDepth;) {
  camPath.Append(sampleBSDF, uniDist);

// Actually sample the path 

Compute an Estimate

We have a path will length maxDepth. In order to create an estimate, we want to take each prefix of our path and connect it to an emitter so that we end up with a complete 3-vertex path, 4-vertex path, etc.

Spectrum Li(0);
for (int length = 2; length < camPath.Size(); length++) {
  // Obtain a prefix subpath
  auto subPath = camPath.Slice(0, length);

  // Append an emitter vertex 
  subPath.Append(sampleEmtPosition, uniDist, scene->emitters);

  // Actually sample the emitter vertex (all the previous vertices have 
  // already been sampled; they will not be resampled)

  // Update estimate
  Li += integrand(subPath) / subPath.Pdf();
return Li;

This loop introduces the Slice method, which can be used on a RandomSequence to obtain a subsequence. The returned subsequence is a valid RandomSequence so it supports all the same methods. Aether also provides methods for concatenating sequences and reversing sequences, which are useful in more complex algorithms.

This example illustrates clearly some of the advantages of Aether: strategies can be easily reused, as can the integrand function; and the RandomSequence data structure makes it easy to assemble paths. The amount of code required to extend direct lighting to full path tracing is relatively small.

Compute an Estimate With MIS

Now we extend the path tracer to use MIS and combine 2 paths: 1 sampled by the BSDF and the other sampled by the light. We only need to make a few changes to the loop body:

// Obtain a subpath sampled by the BSDF
auto bsdfPath = camPath.Slice(0, length + 1);

// BSDF sampled path + sample point on light
auto lightPath = camPath.Slice(0, length);
// Sample the light
lightPath.Append(sampleEmtPosition, uniDist, scene->emitters);

// Combine BSDF path and light path
// Returns a list of paths with their MIS weights
auto combinedList = combine<PowerHeuristic>(bsdfPath, lightPath);

// Sum up the contributions
for (const auto &combined : combinedList) {
  const auto &path = combined.sequence;
  Li += combined.weight * (integrand(path) / path.Pdf());

We obtain 2 subpaths: 1 from sampling the BSDF (which is just camPath) and the other from camPath with 1 fewer vertex, which is then extended by sampling a light vertex. Aether provides the combine primitive to make it easy to compute MIS weights. The PowerHeuristic is provided by Aether. But the user can easily implement alternative combining heuristics. The result of combine is a list of weight, path pairs. Updating the estimate simply involves multiplying the given weight by the standard estimate.