Stratified Sampling

Stratified Sampling

January, 2012

Introduction

RenderMan Pro Server 16.4 introduces the new randomstrat() function to RSL. This new feature simplifies the generation of random samples in one and two dimensions and can reduce image noise when used in conjunction with distribution ray tracing, such as the Physically Plausible shading framework provides. This application note provides details on how to use this new shadeop.


Jittering

Jittering is a widely used stratification technique for reducing noise variance in distribution ray tracing. Typically, we do this by dividing the [0,1) x [0,1) unit square into equally sized cells. For each cell we choose a uniform random point from within it and then warp the point from the unit square to the required domain to be sampled over (e.g., points on the surface of a light). By dividing the unit square into strata and sampling those rather than just uniformly sampling the unit square directly, the amount of sampling noise produced diminishes much more rapidly as the number of samples increases.

Common applications of this include sampling over the surfaces of area lights and sampling over hemispheres to compute BRDFs, ambient occlusion, indirect lighting, and environment lighting. For example, here is a fairly basic rectangular area light source shader:

light arealight( float samples0 = 144;
                 float intensity = 1.0;
                 color lightcolor = 1.0;
                 float falloff = 1; )
{
    normal nn = normalize(N);
    uniform float depth = 0;
    rayinfo("depth", depth);
    uniform float samples = 1;
    if (depth == 0)
        samples = samples0;
    uniform float xstrata = max(floor(sqrt(samples)), 1);
    uniform float ystrata = floor(samples / xstrata);
    uniform float invxstrata = 1.0 / xstrata;
    uniform float invystrata = 1.0 / ystrata;
    samples = xstrata * ystrata;
    point tos[samples];
    color trans[samples];
    vector ls[samples];
    float dists[samples];
    illuminate (Ps + nn) {
        uniform float sample = 0;
        uniform float xstep, ystep;
        for (ystep = 0; ystep < ystrata; ystep += 1)
            for (xstep = 0; xstep < xstrata; xstep += 1) {
                float x = (xstep + random()) * invxstrata;
                float y = (ystep + random()) * invystrata;
                tos[sample] = point "shader" (x - 0.5, y - 0.5, 0.0);
                sample += 1;
            }
        areashadow("", Ps, tos, trans, ls, dists, "raytrace", 1);
        for (sample = 0; sample < samples; sample += 1)
            Cl += intensity * lightcolor *
                  max(ls[sample] . nn, 0) *
                  pow(dists[sample], -falloff) *
                  (1 - trans[sample]);
        Cl /= samples;
    }
}

And here is an example image that uses it:

images/figures.stratification/jitter.png

The two mirrors and the teapot use a surface shader that contains a similar pair of stratification loops for sampling over the hemisphere at the surface and computing an isotropic version of the Ward BRDF.

Note that the two shaders only use a significant number of ray samples when shading points directly visible by the camera. After the first reflection they use only a single sample at each ray hit. Reducing the number of rays at each bounce avoids an exponential increase in the number of rays and keeps render times down. It is quite common to drop down to a single sample or fewer after just a bounce or two. (See The Importance of Importance for more on this topic.)

Unfortunately, when reduced to a single sample per bounce we lose all the benefit of the jitter and are essentially shooting rays using unstratified samples. This is particularly visible in the noisiness of the reflection of the penumbra in the left mirror, compared to its directly visible counterpart. Unless we can somehow stratify across related samples in the ray tree we would need to take either significantly more samples at the first bounce or multiple samples at the later bounces in order to resolve the noise. The new randomstrat() shadeop helps to address this need.


The randomstrat() Shadeop

Given an array and a number of samples to generate, randomstrat() fills in all or part of the array with stratified random samples:

void randomstrat( uniform float sampleoffset, uniform float nsamples,
                  output varying float[] samples );
void randomstrat( uniform float sampleoffset, uniform float nsamples,
                  output varying triple[] samples );
void randomstrat( uniform float sampleoffset, uniform float nsamples,
                  output varying __radiancesample[] samples );

If the array is of type float it produces one-dimensional samples in [0, 1). Otherwise, with an array of points, vectors, normals, or colors it produces two-dimensional samples in the unit square (and zeros out z). It can also fill in the direction field, L, in an array of __radiancesample structs. It fills in the array beginning at sampleoffset and goes until it has produced nsamples. Variable length arrays will be automatically resized if necessary.

Here is the area light shader source from above, modified to use randomstrat():

light arealight( float samples0 = 144;
                 float intensity = 1.0;
                 color lightcolor = 1.0;
                 float falloff = 1; )
{
    normal nn = normalize(N);
    uniform float depth = 0;
    rayinfo("depth", depth);
    uniform float samples = 1;
    if (depth == 0)
        samples = samples0;
    point tos[samples];
    color trans[samples];
    vector ls[samples];
    float dists[samples];
    illuminate (Ps + nn) {
        randomstrat(0, samples, tos);
        uniform float sample;
        for (sample = 0; sample < samples; sample += 1) {
            float x = tos[sample][0];
            float y = tos[sample][1];
            tos[sample] = point "shader" (x - 0.5, y - 0.5, 0.0);
        }
        areashadow("", Ps, tos, trans, ls, dists, "raytrace", 1);
        for (sample = 0; sample < samples; sample += 1)
            Cl += intensity * lightcolor *
                  max(ls[sample] . nn, 0) *
                  pow(dists[sample], -falloff) *
                  (1 - trans[sample]);
        Cl /= samples;
    }
}

And here is the corresponding image (with an isotropic Ward shader similarly modified):

images/figures.stratification/randomstrat.png

Though the ray counts and render times are equivalent, the noise in the reflections has been reduced, particularly in the sharp reflection on the left. The following image shows detail of a couple of these areas. For each pair, the RSL jittered version is on the left and the randomstrat() version is on the right:

images/figures.stratification/comparison.png

Under the Hood

Internally, PRMan can now keep track of which rays are related by common ancestry in the ray tree. As the rays propagate from bounce to bounce they inherit this information. The randomstrat() shadeop can use this information to provide repeatable samples that are stratified across multiple rays.

Here is how PRMan currently decides which rays are related. If you invoke a ray tracing shadeop such that it shoots multiple rays per shading point, we consider the rays in each of those sets to be related amongst themselves and independent from the other sets of rays. For example, if we shoot five rays on the first bounce and then three rays each on the second bounce, the ray tree looks like this where each unique color shows rays considered to be related for stratification purposes:

images/figures.stratification/broad.png

On the other hand, if a ray tracing shadeop is used to shoot only a single ray at a time, such as in path tracing, we inherit the grouping from the parent rays. Additional bounces that extend the ray paths will continue to inherit the grouping from the parent rays:

images/figures.stratification/narrow.png

Suppose that our shader had fired two rays, each from a separate call to a ray tracing shadeop. In this case, all the rays from the first call would be considered to be related and any single-ray bounces that extend those paths would also be seen as related. Any calls to randomstrat() or the ray tracing system after that would be considered to create a new group of rays. Those rays and their decendents would still be stratified between themselves, but the new paths would be fully independent of the first paths:

images/figures.stratification/split.png

In short, sending a batch of rays from a single shading point always creates a new stratification group, regardless of any previously inherited grouping. Sending a single ray per shading point causes a new group of paths that inherit the stratification of their parent rays.

The randomstrat() shadeop itself works similarly. When nsamples is greater than one, it ignores the existing stratification groups and generates a new batch of samples. These will be internally stratified but the set at each shading point will be independent of the sets at the other shading points. When nsamples is one, the sample returned in the ray will be stratified across each related set of shading points as determined by the ray tree inheritance.

Note that in general, the ray tracing shadeops determine the stratification relationship that randomstrat() will assume. For best results with randomstrat(), you will typically want to have a one-to-one correspondence between an invocation of it and a subsequent launching of rays.


Future Directions

  • Only explicit sample generation in RSL via randomstrat() produces stratification across related rays when tracing a single sample at a time from each shading point. The generateSamplesAS(), generateSamplesEnv(), gather(), occlusion(), and similar shadeops do not yet take advantage of this new capability.
  • In the current release, we distinguish between stratification over single vs. multiple samples at a time when invoking randomstrat() or ray tracing. In the future, we would like to allow cooperative stratification. For example, in the case above of a shading point tracing five rays on the first bounce, each of which traces three on the second bounce, we would like to have each group of three be stratified but in a such a way that all fifteen rays are also well stratified.