Domain Sampling

OpenQMC aims to provide an API that makes using high quality samples easy and bias-free when writing software with numerical integration. You have already seen a very basic example of this in the Usage section. This section will go into more detail on the concepts and provide some more complex examples.

The section will build upon a simple example, introduce the concepts of domain trees and sample splitting, and show some of the potential pitfalls OpenQMC prevents by following some basic rules. After this, you should be able to confidently make use of QMC sampling in your own software.

Sampling domains

When you estimate an integral using the Monte Carlo method, the high dimensional space can often be partitioned into logical domains. In OpenQMC domains are an important concept. Domains act as promises that can be converted into dimensions upon request by the caller. They can also be used to derive other domains.

This example function tries to compute a pixel estimate. This estimate is made up of two parts, each needing a domain. First is a camera that takes a single dimension, and the second is a light, that takes two dimensions.

Result estimatePixel(const oqmc::PmjBnSampler pixelDomain,
                     const CameraInterface& camera,
                     const LightInterface& light)
{
    enum DomainKey
    {
        Next,
    };

    // Derive 'cameraDomain' from 'pixelDomain' parameter.
    const auto cameraDomain = pixelDomain.newDomain(DomainKey::Next);

    // Derive 'lightDomain' from 'cameraDomain' variable.
    const auto lightDomain = cameraDomain.newDomain(DomainKey::Next);

    // Take a single dimension for time to pass to 'camera'.
    float timeSample[1];
    cameraDomain.drawSample<1>(timeSample);

    // Take two dimensions for direction to pass to 'light'.
    float directionSample[2];
    lightDomain.drawSample<2>(directionSample);

    // Pass drawn dimensions to the interface for importance sampling.
    const auto cameraSample = camera.sample(timeSample);
    const auto lightSample = light.sample(directionSample);

    return estimateLightTransport(cameraSample, lightSample);
}

The function started with an initial pixelDomain, from that it derived a cameraDomain, and from that it derived a lightDomain. This is a linear sequence of dependent domains. Such an operation is called chaining.

Domains provide up to 4 dimensions each, often referred to here as a sample (technically part of a sample). If more than 4 dimensions are required, the caller can always derive a new domain.

This example works, but isn't very flexible. This can be improved. Types that implement the camera or light interfaces may require less or more dimensions than is provided by the calling code.

Passing domains

Computing domains is cheap, but drawing samples is expensive. In the last example, a sample was drawn and passed to a type that implements the LightInterface. If the sample is not required, the type doesn't have the option to prevent the sample from being drawn.

You can resolve this by changing the interface to pass a domain, deferring the decision to draw a sample to the type. Here are two types that implement such an interface, DiskLight and PointLight.

LightSample DiskLight::sample(const oqmc::PmjBnSampler lightDomain)
{
    // Draw two dimensions for direction.
    float directionSample[2];
    lightDomain.drawSample<2>(directionSample);

    // Compute LightSample object using the drawn sample...
}

LightSample PointLight::sample(const oqmc::PmjBnSampler lightDomain)
{
    // Don't draw two dimensions for direction.
    // float directionSample[2];
    // lightDomain.drawSample<2>(directionSample);

    // Compute LightSample object without using a drawn sample...
}

The DiskLight draws a sample from the domain. But PointLight opts to save on compute and does not draw a sample as it is not required.

Looking at the revised calling code, the domains are now passed through the interfaces, allowing the types to decide whether to draw a sample or not.

Result estimatePixel(const oqmc::PmjBnSampler pixelDomain,
                     const CameraInterface& camera,
                     const LightInterface& light)
{
    enum DomainKey
    {
        Next,
    };

    // Derive 'cameraDomain' from 'pixelDomain' parameter.
    const auto cameraDomain = pixelDomain.newDomain(DomainKey::Next);

    // Derive 'lightDomain' from 'cameraDomain' variable.
    const auto lightDomain = cameraDomain.newDomain(DomainKey::Next);

    // Pass domains directly to the interface to compute a sample.
    const auto cameraSample = camera.sample(cameraDomain);
    const auto lightSample = light.sample(lightDomain);

    return estimateLightTransport(cameraSample, lightSample);
}

Sometimes a type needs an unknown number of domains. The next section will show how that can be achieved. This will introduce a potential danger, as well as a new concept called domain trees that can be used to do the operation safely.

Domain branching

Building on the previous example, the calling code is free of knowing how many domains a type might need. If a type requires more sample draws, it can always derive a new domain from the domain that was passed to it, and then draw more samples from those domains.

This is an example where ThinLensCamera implements the CameraInterface. It requires multiple domains and uses chaining to sequentially derive domains one after the other, starting with the original cameraDomain.

CameraSample ThinLensCamera::sample(const oqmc::PmjBnSampler cameraDomain)
{
    enum DomainKey
    {
        Next,
    };

    // Derive 'lensDomain' from 'cameraDomain' parameter.
    const auto lensDomain = cameraDomain.newDomain(DomainKey::Next);

    // Derive 'timeDomain' from 'lensDomain' variable.
    const auto timeDomain = lensDomain.newDomain(DomainKey::Next);

    // Take two dimensions for lens to compute a CameraSample object.
    float lensSample[2];
    lensDomain.drawSample<2>(lensSample);

    // Take a single dimension for time to compute a CameraSample object.
    float timeSample[1];
    timeDomain.drawSample<1>(timeSample);

    // Compute CameraSample object using the drawn samples...
}

Here is the dangerous part! Notice how the calling estimatePixel function derives both lensDomain and lightDomain from cameraDomain, meaning they become equivalent in value.

This will result in lens and light dimensions correlating, leaving gaps in the primary sample space. That type of correlation will bias your estimates, and should be avoided. A revised version of the estimatePixel function:

Result estimatePixel(const oqmc::PmjBnSampler pixelDomain,
                     const CameraInterface& camera,
                     const LightInterface& light)
{
    enum DomainKey
    {
        Camera,
        Light,
    };

    // Derive 'cameraDomain' from 'pixelDomain' parameter.
    const auto cameraDomain = pixelDomain.newDomain(DomainKey::Camera);

    // Derive 'lightDomain' from 'pixelDomain' parameter.
    const auto lightDomain = pixelDomain.newDomain(DomainKey::Light);

    // Pass domains directly to the interface to compute samples.
    const auto cameraSample = camera.sample(cameraDomain);
    const auto lightSample = light.sample(lightDomain);

    return estimateLightTransport(cameraSample, lightSample);
}

You may notice that the DomainKey enum has now changed. The new code derives both cameraDomain and lightDomain directly from pixelDomain, but in each instance passes a different enum value.

This operation is called branching. The derived domains are now branched and independent. You can branch a domain on a given key either with an enum as seen here, or by passing an integer value directly.

Now that cameraDomain and lightDomain are independent, the ThinLensCamera can use the domain passed to it in whichever way it needs.

Mapping domains to the code's call graph creates domain trees. If you carefully derive domains by using chaining when writing loops, and branching when passing to interfaces, you'll avoid gaps in the primary sample space, guaranteeing bias-free results.

OpenQMC's design makes domain trees a first-class concept and lets you focus on writing complex and correct sampling code by applying the techniques described.

Here is what this example has identified. Gaps in the primary sample space produce bias. Branching can be used to make domains independent of one another. This independence prevents such bias. Finally, domain trees should match the call graph of the code to guarantee bias-free results. For a more complete example, see the trace tool.