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.

MIS choices in the first HdRestir commit

PBRT §2.2.3 supplies the general MIS derivation. The relevant local question is how the initial commit encodes its density assumptions and numerical guards.

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.

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 expectation-preserving rescaling 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

Rendered convergence: one fixed random sequence Move the control to inspect the same Path Tracer scene from 1 to 64 samples per pixel.

Samples per pixel

8 spp

8 spp

The displayed run becomes less noisy overall, but individual pixels need not approach their limiting values monotonically. Moving from \(N\) to \(4N\) samples halves the standard-deviation scale under the independence assumptions derived above; it does not divide every visible error by two.

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