An embedded domain specific sampling language for Monte Carlo rendering

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

Aether is available here.

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:

```
-stdlib=libc++
-std=c++1z
-ftemplate-depth=512
-fconstexpr-steps=127124200
-fconstexpr-depth=720
```

All Aether code is wrapped in the `aether`

namespace.

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

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.

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]);
```

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.

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 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 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);
```

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.

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
path.Sample();
```

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.

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.

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.

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;
```

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.

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(
pt
, 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);
```

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(
pt
, 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);
```

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(
pt
, 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();
trianglePtrs.reserve(triangleCount);
for (size_t i = 0; i < triangleCount; i++) {
trianglePtrs.push_back(trianglePtr++);
}
triangleWeights.reserve(triangleCount);
for (size_t i = 0; i < triangleCount; i++) {
// Sample according to area
triangleWeights.push_back(Real(trianglePtrs[i]->surfaceArea(vertexPositions)));
}
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.

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:

```
camPath.Sample();
```

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`

:

- We can access the _i_th sample with
`path[i]`

- The actual 3D vertex of the _i_th sample can be accessed with
`path[i].Value()`

- The number of samples in the path is given by
`path.Size()`

- Additional data for the _i_th sample can be accessed with
`path[i].Get(intersection_)`

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:

```
camPath.Pdf();
```

So the final estimate is simply:

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

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);
```

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 ¤tIts = 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(
nextIts.isValid()
, 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
camPath.Sample();
```

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)
subPath.Sample();
// 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.

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);
lightPath.Sample();
// 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.