# 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:

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:

voidrandomstrat( uniform float sampleoffset, uniform float nsamples, output varyingfloat[]samples ); voidrandomstrat( uniform float sampleoffset, uniform float nsamples, output varyingtriple[]samples ); voidrandomstrat( 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):

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:

### 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:

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:

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:

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.