Monte Carlo path tracing

Path tracing evaluates the rendering equation with Monte Carlo integration. The ray sequence is an implementation of a recursive integral: each intersection supplies the next integration domain, and each sampled direction contributes a factor consisting of a BSDF value, a cosine, and a reciprocal probability density.

This chapter extends the site’s earlier Monte Carlo introduction from volumes to surfaces.

The rendering equation

CUDA VolPT introduces radiance and the radiative-transport problem for participating media. The equation below is the surface-scattering counterpart: both describe how emitted and transported radiance combine, but the integration domain here is the hemisphere of directions at a surface.

At surface point \(x\), outgoing direction \(\omega_o\), and wavelength \(\lambda\),

• \(L_e\): light emitted by the surface.
• \(L_i\): radiance arriving from direction \(\omega_i\).
• \(f_s\): the BSDF, describing how the material redirects light.
• \(\lvert n\cdot\omega_i\rvert\): projected-area cosine.
• \(\Omega^+\): the hemisphere above the surface.

The equation is recursive because arriving radiance is outgoing radiance from another point. A path tracer samples a direction, traces to that point, and repeats until the path terminates.

The terms must use compatible measures. \(f_s\) has units of inverse steradians, \(d\omega_i\) has units of steradians, and the cosine is dimensionless. Their product therefore leaves the units of radiance unchanged. This unit check is useful when reviewing code: adding an area-density PDF to a solid-angle-density PDF is not only mathematically incorrect; their units are different.

For an opaque surface the hemisphere is normally oriented by the shading normal. Transmission introduces directions on the opposite side, and delta distributions represent ideal reflection or refraction. HdRestir’s material sampling code handles those cases through concrete BSDF implementations and flags; the integral above is the non-delta form used to establish the basic estimator.

The Monte Carlo estimator

Write the integral as \(I=\int_D f(x)\,dx\). Draw \(X_i\sim p(x)\):

The variance statement describes the distribution of repeated estimates. It does not say that the absolute error of one random sequence decreases after every additional sample. The following model estimates a one-dimensional integral with a known exact value so those two claims can be separated.

The division by \(p\) corrects for sampling frequency. Sampling a region more frequently without dividing by its density would change the integral being estimated. Conversely, a very small \(p(x)\) where \(\lvert f(x)\rvert\) is large produces large weights and high variance even though the estimator remains unbiased.

Importance sampling

Variance disappears completely in the ideal positive case \(p^\star(x)=f(x)/I\), because \(f(X)/p^\star(X)=I\) for every sample. We do not know \(I\), and sampling the full lighting integrand is usually hard, but we can approximate its shape:

• BSDF sampling follows directions favored by the material.
• Light sampling chooses points or directions on emitters.
• Environment sampling follows bright regions of an HDR map.

ReSTIR does not remove this requirement. It constructs a discrete resampling distribution over generated candidates and applies a second weight that corrects for that selection. Chapters 6 and 9 derive that correction.

Direct lighting and visibility

Next-event estimation (NEE) explicitly samples a light at a surface. For a finite area light sampled by area, a point \(y\) with area density \(p_A(y)\) must be converted to solid-angle density at \(x\):

The factor is the area-to-solid-angle Jacobian. It is not optional: the BSDF and the rendering equation are measured in solid angle. A shadow ray supplies the binary visibility \(V(x,y)\).

The first commit already performs next-event estimation, but it stores the light and BSDF densities as unlabelled float values:

const float totalLightPdf{lightSelectPdf * ls.Pdf};
const float misWeight{
    light.IsDeltaLight()
        ? 1.0f
        : PowerHeuristic(totalLightPdf, bsdfPdf)};

This literal excerpt comes from direct_lighting.cpp at the first commit. lightSelectPdf is the probability mass for selecting a light, while ls.Pdf is the conditional density returned by that light. Their product is the joint proposal density. The code assumes that ls.Pdf and bsdfPdf are already comparable for a non-delta light; the type system does not verify the measure.

The explicit PdfSpace type and area-to-solid-angle conversion are introduced by the next commit and are examined in step 0.g. They are not projected backward into this foundation chapter because they do not exist in the first revision.

If shadow maps are new to you, the site’s shadows tutorial builds intuition for visibility; a path tracer answers the same question with a ray instead of a depth-map lookup.

Multiple importance sampling

A glossy BSDF is good at finding a small reflection lobe; light sampling is good at finding a small emitter. Multiple importance sampling (MIS) combines both. Suppose techniques \(t=1,\ldots,T\) draw \(N_t\) samples with densities \(p_t\). The balance heuristic uses

and the estimator

This is the central result of Veach’s MIS formulation. HdRestir uses MIS between light, environment, and BSDF proposals.

Power heuristic in the implementation

The first commit implements the two-technique power heuristic directly:

inline float PowerHeuristic(float f, float g) {
    float f2 = f * f;
    float g2 = g * g;
    return f2 / (f2 + g2);
}

For one sample from each technique this is \(p_f^2/(p_f^2+p_g^2)\). The helper has a numerical precondition not encoded by its signature: at least one argument must be positive and finite. If both are zero, the denominator is zero. The direct-light caller avoids applying the continuous heuristic to delta lights, but every continuous call site must still satisfy this precondition.

Direct-light contribution in the first commit

The rest of the first-commit expression is:

return RGBToSpectrum(bsdfVal, surface.lambda)
     * RGBToSpectrum(ls.Color, surface.lambda)
     * (nDotL / (totalLightPdf + 1e-6f))
     * misWeight;

The BSDF, emitted light, and cosine form the numerator. The joint light PDF is in the denominator, and the MIS weight partitions the result with the BSDF-sampled technique. The added \(10^{-6}\) is a numerical regularizer, not part of the exact estimator. Replacing \(p\) with \(p+\varepsilon\) changes the expectation for arbitrary \(p\). Later chapters apply the same distinction between an exact derivation and a historically convenient denominator guard.

Paths and Russian roulette

The thesis’s volumetric path-tracer algorithm shows the same throughput structure with an additional distance-sampling step inside participating media. HdRestir uses the surface branch of that construction at every mesh intersection.

At bounce \(k\), the path throughput multiplies by

Products of these ratios carry the contribution of the whole path. To avoid tracing negligible paths forever, Russian roulette continues with probability \(q\). If continued, throughput is divided by \(q\).

Path construction in the first commit

There is no PathIntegrator class yet. The bounce loop is a local TracePath function in path_trace_pass.cpp. The following excerpt is literal; unrelated opacity and absorption branches are marked as omitted:

SampledSpectrum throughput{1.0f};
SampledSpectrum totalRadiance{0.0f};
Ray currentRay{cameraRay.origin, cameraRay.dir};
std::optional<HitRecord> nextHit{primaryHit};
BounceState bounceState{};

for (int bounce{0}; bounce < maxDepth; ++bounce) {
    if (!nextHit.has_value()) {
        if (env != nullptr && (bounce > 0 || settings.RenderIblBackground)) {
            totalRadiance += throughput * RGBToSpectrum(env->Sample(currentRay.Dir), lambda);
        }
        break;
    }

    HitRecord hit{*nextHit};
    const IMaterial* material{scene.GetMaterial(hit.MatId)};
    if (material == nullptr) {
        material = &DefaultMaterial::Instance();
    }

    // ...

    totalRadiance += throughput * RGBToSpectrum(c.Emission, lambda);
    const std::unique_ptr<IBSDF> bsdfOwner{material->CreateBSDF(BSDFClosure{c})};
    const ShadingPoint surface{hit, *bsdfOwner, c, shadingNormal, currentRay.Dir, lambda, isInside};
    totalRadiance += throughput * SampleDirectLighting(surface, scene.GetLights(), scene, rng);

    const BounceConfig config{settings.MaxReflectionBounces, settings.MaxRefractionBounces};
    const BounceSample bs{material->SampleBounce(surface, config, bounceState, rng)};

    if (bs.Terminate) {
        break;
    }

    throughput *= bs.ThroughputMul;
    currentRay = bs.NextRay;
    nextHit = scene.IntersectScene(currentRay.Origin, currentRay.Dir);

    // ...
}

throughput is the product of all preceding bounce factors. Emission and direct light at the current surface are multiplied by that product before a new bounce is sampled. BounceSample::ThroughputMul already contains the BSDF, cosine, and reciprocal sampling density for the selected event. The refactor that moves this loop into PathIntegrator belongs to commit e6fb9bb and is therefore left to step 0.g.

Exact Russian-roulette condition

The first commit contains this code:

if (bounce > 3) {
    const float p{throughput.Max()};
    if (rng.NextFloat() > p) {
        break;
    }
    throughput *= 1.0f / p;
}

It expresses the intended correction but does not enforce the mathematical domain of a probability. throughput.Max() can be zero or greater than one. Moreover, NextFloat() can return zero. If p == 0 and the draw is zero, the comparison is false and the following division is undefined. If p > 1, the path always survives but is nevertheless divided by a number greater than one, which changes its expectation.

A mathematically valid form is the following derived code, not a quotation from the commit:

const float q{std::clamp(throughput.Max(), 0.0f, 1.0f)};
if (q == 0.0f || rng.NextFloat() >= q) {
    break;
}
throughput *= 1.0f / q;

For \(q=0\), a zero-throughput path terminates before division. For \(0<q<1\), survival has probability \(q\) because a uniform value in \([0,1)\) satisfies \(U<q\) with probability \(q\). For \(q=1\), the path always survives and division leaves throughput unchanged. These three cases satisfy the assumptions stated in the proof for finite, non-negative path throughput. Non-finite values indicate an earlier numerical failure and require separate rejection or diagnostics; they are not valid roulette probabilities.

Implementation files

The next chapter explains the Hydra boundary through which usdview supplies the scene and receives the buffers produced by this first implementation.