Resampled Importance Sampling

This commit adds a direct-light integrator based on Resampled Importance Sampling (RIS). It generates several candidates from the existing light and BSDF proposals, computes a scalar weight for each candidate, retains one candidate by weighted categorical sampling, and applies the RIS normalization to its full spectral contribution.

The idea originates in Talbot, Cline, and Egbert’s Importance Resampling for Global Illumination.

The algorithm

Let \(f(x)\) be the contribution we want to integrate. Draw \(M\) candidates \(X_i\sim q(x)\). Choose a positive scalar target \(\hat p(x)\), which should track the size of \(f\), and compute

Select \(Y=X_i\) with conditional probability \(w_i/S\). The RIS estimate is

The target need not be normalized. Its unknown constant cancels.

Scaling the target by any constant \(a>0\) multiplies every \(w_i\) and \(S\) by \(a\), leaves the categorical probabilities \(w_i/S\) unchanged, and cancels in \(S/(M\hat p(Y))\). The target controls relative selection within a candidate set; its absolute normalization is irrelevant in exact arithmetic.

The model below uses four fixed candidates. The target-scale control multiplies every unnormalized target by the same constant. The weight sum changes, while the selection probabilities and final correction do not. The selection control chooses a point on the resulting categorical distribution.

For \(M=1\), RIS reduces to ordinary importance sampling:

This is a useful implementation test. Running the RIS pipeline with one candidate should match the corresponding proposal estimator in expectation. It need not produce the same individual image if the two paths consume random numbers in a different order.

flowchart LR
  Q["proposal q"] --> C["M candidates"]
  C --> W["weight p̂/q"]
  W --> Y["retain one<br/>probability ∝ weight"]
  Y --> E["evaluate / shade one sample"]
  W --> E

What the target should contain

The ideal target is proportional to contribution magnitude, but it must be cheap. HdRestir scalarizes its four-component spectral contribution with luminance. Candidate weights also include the relevant MIS factor and proposal PDF. The selected sample keeps spectral emission and throughput; the scalar target decides which candidate survives, not the final color.

This distinction prevents a luminance target from turning the transport itself into monochrome.

The target affects variance. If it is nearly proportional to the magnitude of the final contribution, high-value candidates are retained frequently and the correction varies less. If it predicts contribution poorly, the algorithm can select samples with a large target but small final value, or rarely select a large contribution and assign it a large correction. The estimator can remain unbiased in both cases while their convergence rates differ substantially.

Multiple proposal techniques

HdRestir generates NEE and BSDF candidates. Each candidate’s proposal density and MIS factor must refer to the technique that produced it. One safe view is to treat “technique index + sample” as the sampling domain. The same conditional proof applies after replacing \(q\) by the correct joint proposal.

Let \(T\) be a discrete technique index with selection probability \(\pi_t\), and let \(q_t(x)\) be the conditional proposal. The joint proposal is

Every RIS weight must divide by this joint density unless the technique selection and sample-count factors have already been incorporated through the MIS coefficient. This is the reason the code review must trace all normalization factors across candidate construction rather than inspecting only the final return statement.

The exploratory programs ris_1d.py and ris_mis_two_techniques.py test known integrals before the estimator is trusted in 3D.

The running mean approaching the analytic result checks bias empirically. It does not replace the proof; it catches implementation errors that the proof’s assumptions do not.

Read the RIS diff in estimator order

The implementation matches the derivation above almost line for line: construct candidates, assign scalar weights, sample a categorical index, then apply the correction that makes the retained contribution unbiased.

NEE candidate weight

const MISContrib nee{
    _evaluateNEE(surface, *candidate->Light, candidate->Ls, scene)};
const float targetFunction{
    SpectrumLuminance(
        nee.Emission * nee.ThroughputIntegrandMul, surface.lambda)};
const float misWeight{(nee.UseMis && useBsdfTechnique)
    ? PowerHeuristic(
          nee.PNee, candidateCount, nee.PBsdf, candidateCount)
    : 1.0f};

// ...

const double risWeight{
    misWeight * targetFunction / (nee.PNee + 1e-9)};

Read the last line as \(w=m\,\hat p/q\):

• targetFunction is \(\hat p\), a scalar luminance proxy;
• nee.PNee is proposal density \(q\);
• misWeight partitions contribution between NEE and BSDF techniques.

The spectral Emission and ThroughputIntegrandMul remain in the candidate. Only selection uses the scalar target.

The 1e-9 term is a numerical regularizer. It prevents division by zero when a reported PDF is zero, but it is not part of the exact RIS derivation:

Therefore the implementation with nonzero \(\varepsilon\) is not exactly the unbiased mathematical estimator for arbitrary \(q\). The difference is normally negligible when valid candidate PDFs are much larger than \(10^{-9}\), but a stricter implementation would reject candidates whose PDF is non-positive and divide positive PDFs without adding epsilon. Tests should include very small PDFs so this behavior is explicit.

Weighted categorical selection

double remainingWeight{
    static_cast<double>(rng.NextFloat()) * weightSum};
for (int i{0}; i < gsl::narrow<int>(candidates.size()); ++i) {
    const double weight{
        candidates[static_cast<std::size_t>(i)].GetWeight()};
    if (weight <= 0.0) {
        continue;
    }
    remainingWeight -= weight;
    if (remainingWeight <= 0.0) {
        return ChosenCandidate{
            .Index = i,
            .WeightSum = weightSum};
    }
}

A uniform number selects a point on \([0,S)\). Subtracting each positive weight walks through adjacent intervals of lengths \(w_i\); candidate \(i\) therefore owns exactly \(w_i/S\) of the interval. Returning WeightSum is essential because the final correction needs \(S\), not only the chosen index.

The helper ignores non-positive weights. It should also reject non-finite weights: NaN breaks ordered comparisons and infinity can make rng.NextFloat() * weightSum undefined. The later reservoir implementation adds stronger finite-value guards on the GPU. At this revision, candidate generation is responsible for producing finite target values and PDFs.

Spectral contribution and RIS correction

const RISLightCandidate& chosenCandidate{
    chosenCandidateOptional.value()};

// ...

const float unbiasedContributionWeight{
    static_cast<float>(
        chosen->WeightSum / chosenCandidate.targetFunction)};

return chosenCandidate.Emission
     * chosenCandidate.ThroupoutMult
     * unbiasedContributionWeight;

This checkpoint stores per-technique \(1/M\) factors in its candidate/MIS construction, so the displayed correction appears as WeightSum / targetFunction rather than a literal WeightSum / (M * targetFunction). When reading estimator code, always locate where sample-count normalization lives before deciding that a factor is missing.

This sample-count placement is specific to the implementation at this commit. The mathematical reference remains \(S/(M\hat p)\). Refactoring the MIS helper or changing the number of samples per technique requires re-deriving the code coefficient; copying the same final expression into a new convention can remove or duplicate \(1/M\).

RIS pipeline construction

passSpecs.push_back(
    Detail::MakePassSpec<RISPathTracePass>(
        settings.CandidateCount,
        settings.PathTracer.PathTrace,
        settings.PathTracer.MaxDepth));
passSpecs.push_back(
    Detail::MakePassSpec<AccumulationPass>(
        settings.PathTracer.Denoiser.EnableFireflyFilter));

Only the integration pass changes. Accumulation, optional upscale, AOV selection, and Hydra handoff are shared with the reference path tracer as a result of the interface introduced by the previous commit.

RIS is not yet ReSTIR

At this checkpoint the candidate list exists only for the current shading operation. The losing candidates are discarded and nothing crosses a pixel or frame boundary. RIS supplies the resampling estimator; ReSTIR adds compact reservoir state and reuse.

The cost of this implementation grows approximately linearly with candidate count because it stores and evaluates a vector of candidates before selecting one. Increasing \(M\) can improve the retained distribution but also increases BSDF evaluation, light sampling, and—in this revision—visibility work. The correct comparison is error against elapsed time for several candidate counts, not the assumption that a larger set is always preferable.

Code map and reference

• ris_direct_light_integrator.cpp
• random_index.h
• Original ReSTIR paper, especially its streaming RIS derivation